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PREFACE 


In the past decade or so, there has been considerable interest and progress in the 
development and utilization of space techniques for precise measurements of geodetic 
baselines, earth orientation, and various geodynamic studies, especially for measuring 
large-scale distortions within plates and determining the rates of interplate motion. 
These methods rely heavily on extra-terrestrial reference sources such as the distant 
quasars or other compact extragalactic objects used in Very Long Baseline 
Interferometry (VLBI), or the moon used for Lunar Laser Ranging (LLR) and the low- 
earth satellites such as LAGEOS and STARLETTE used in Satellite Laser Ranging (SLR). 
Notably, VLBI and SLR have achieved significant superiority over other conventional 
approaches for measuring vector baselines very precisely. As currently applied, VLBI 
and SLR have reached a level of maturity that to date can be used to measure routinely 
baseline vectors with lengths up to intercontinental distances with repeatabilities of 
0.01 parts per million or better in both length and orientation. 

In both techniques a number of fixed stations are used to determine the variations 
of the earth’s angular orientation in space, to measure plate tectonic motions through 
monitoring of the locations of the fixed stations with respect to each other, and to 
contribute to the maintenance of a reference frame with respect to which the motion of 
additional points of interest can be determined by means of mobile VLBI and SLR 
equipment. However, in addition to the high cost of instrumentation and operation, such 
mobile systems can still be somewhat limited in their ability to occupy sites which are 
not easily accessible, thus limiting their use for many regional geodetic and geodynamic 
applications where more measurements of this type are needed, at more frequent 
intervals in time and space. Operational costs are particularly high for SLR due to the 
system’s susceptibility to weather. Typically, 5 to 30 days for fixed (and up to 60 days 
for mobile) site occupations are required if the length of intersite baselines up to 
intercontinental distances were to be determined to a precision of 3-5 cm. By contrast, 
for the most basic IRIS network of fixed VLBI sites (i.e. the POLARIS and Wettzell 
observatories) the time typically required to achieve sub-decimeter accuracies in the 
determination of baselines of similar lengths corresponds to observation intervals of the 
order of 24 hours. Mobile VLBI systems are less sensitive to adverse weather, but 
involve considerable operations since the smaller-diameter antennas yield less- 
sensitive interferometers than normally achievable with the larger fixed antennas 
which, in turn, impose severe limitations in the observing schedules (often restricting 
observations to the stronger sources) and tend to distort the experimental geometry and 
observing strategies. This decrease in sensitivity can, in principle, be compensated for 
by high-gain antennas, low-noise receivers and multi-observing sessions but not 
without the expense of all the attendant complexities and increased cost of operations. 

Although these technologies are becoming an increasingly important tool for 
geodynamic studies, the future role of mobile VLBI and SLR may well be fulfilled by 
using alternative techniques such as those utilizing the signals from the Global 
Positioning System (GPS) which, already without the full implementation of the system, 
offers a favorable combination of cost and accuracy and has consistently demonstrated the 
capability to provide high-precision densification control in the regional and local areas 
of the VLBI and SLR networks. Although GPS itself is still technically in its testing phase 
and is not expected to become fully operational until the early 1990’s with the 



placement of 18 operational and 3 active spare satellites in equally-spaced orbits (three 
satellites in each of six orbital planes), it has already proved its usefulness in 
measuring relative and absolute positions time and again. Numerous studies, tests, 
comparisons and actual full-scale projects have shown that currently accuracies of a 

few parts in 10 6 over distances up to a few hundred kilometers and observing time 
periods ranging from of 3-5 hours to as little as a few minutes are being achieved 

routinely with standard receiver equipment. Accuracies of a few parts in 10 7 over 
longer distances are becoming increasingly likely with extended modelling of the 
dominant error sources, whereas with carefully designed experiments, these errors 
may be reduced with special efforts to improve the satellite ephemerides to a few parts 

in 10 8 . This level of performance requires, in addition to the special orbit refinement 
efforts, high-performance GPS receivers and reliable atmospheric calibrations 
(especially of the wet troposphere, for instance, by water-vapour radiometers). 

The high-performance, low-cost versatility and potential of GPS-based geodetic 
systems have sparked an intense interest within several government agencies in the 
U.S., Canada, Europe and Australia which have initiated several system studies and 
identified promising geodetic applications, including the expected modes of operation. 
Performing simultaneous observations at widely separated permanent control sites 
together with mobile GPS receivers would conceivably allow centimeter-level 
accuracies for baselines up to thousands of kilometers in length. Even back-packable 
GPS satellite observing systems would thus afford the ability to measure baselines 
between VLBI or SLR fixed sites and arbitrary locations accessible to the GPS systems. 
This combination will allow several more points to be tied to the VLBI and SLR networks 
with comparable accuracy for detailed regional and local monitoring that is rarely 
available today. This premise has led the U.S. National Aeronautics and Space 
Administration (NASA) and the Jet Propulsion Laboratory (JPL) to explore the 
feasibility of the fiducial concept approach as it is commonly referred to (e.g. Davidson 
et al, 1985) which, as envisaged, will enable the joint determination of GPS satellite 
orbits and geodetic baselines in a manner much more suitable for the kind of 
measurements needed for continuous tectonic monitoring of areas of geophysical and 
geodynamic interest. Similar plans currently underway in Canada (e.g. Delikaraoglou et 
al, 1986) include the implementation of an Active Control System (ACS) based on the 
concept of a sparse network of control stations where continuous GPS tracking is carried 
out. This GPS network would be combined with a number of widely spaced, fixed VLBI 
antennas to ensure consistent scale and orientation of the terrestrial networks positioned 
with GPS with respect to the ACS system. 

There are several factors that must be considered in the assessment of the 
potential applications and the need for a strategy towards the implementation of a fully 
operational fiducial or ACS network. There seem to be three main tasks in the way to the 
full implementation of such a system: setting up automated GPS stations (e.g. the 
Canadian ACS units are designed to be microcomputer-controlled, with suitable 
communication interfaces and ability to monitor the station operations and accomplish 
other station-keeping functions automatically); financing the whole operation; and 
analyzing the data to generate and distribute the satellite orbits and organize an efficient 
system and means of transmitting ACS data to the various users. All that can be stated 
about the first two is that they are quite considerable. The last one is also on the 
formidable side; of the three, this is the only one which this work is concerned with. 
The main approach has been based on three main issues, that firstly, to obtain accuracies 
at 0.01 ppm or better, existing models and software for the processing of GPS 
differential observations should be improved considerably through analysis of existing 
and future GPS observations; furthermore, investigations in improving the accuracy of 
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the orbital information for the GPS satellites should continue through modelling 
refinements of the orbital dynamics of this particular type of satellite; and finally, 
consideration of the problem of the optimal GPS network design strategies should be 
undertaken. 

In assembling the material for this report we have reviewed in Section 1 the need 
for both VLBI and SLR vis-&-vis GPS and have outlined the capabilities and limitations of 
each technique and how their complementary applications can be of benefit to geodetic 
and geodynamic operations. To establish the accuracy requirements using GPS for a 
relative positioning accuracy at the level of 0.01 ppm or better, the effects of various 
biases and errors are briefly examined in Section 2. Ways of minimizing these effects 
in general and within the context of this study in particular are also examined. Of these, 
the orbit errors have been recognized since the early days of GPS as constituting a 
severe problem. Things have improved considerably lately, although these errors still 
remain mostly a nuisance that can be mitigated with various sophisticated models, but 
still limit the usefulness of the data in many applications. For this reason, the 
development of the models pertinent to the problem of orbit estimation for the GPS 
satellites receives particular attention in Section 3, where some detailed background and 
the models used in this study are presented. In doing so, we have followed the standard 
formalism based on celestial mechanics which can give further insight into the way 
empirical methods work and perhaps into how to improve them, in the hope that this 
rather detailed introduction to the otherwise basic concepts can be of some value to those 
who wish to understand the nature of the orbit errors better and how these are shaped 
for the particular type of GPS satellite. The basic principles of the coordinate and time 
frames inherent in the reduction of GPS data are reviewed in Section 4. Current 
activities in establishing GPS-based automated geodetic control systems are reviewed in 
Section 5. Strategies for the simultaneous determination of GPS satellite orbits and 
geodetic baselines are also examined in Section 5, followed by actual results of data 
reductions carried out to test these strategies and the general methodology described in 
this report. Finally, a summary of the goals achieved with this study and directions for 
future work are given in Section 6. 
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SECTION 1 


THE NEED FOR ULBI/SLR UIS-R-UIS GPS 


1.1 Requirements for Geodetic and Geodtjnamic Control 


Geodesy and geodynamics derive much of their observational techniques from 
precise positioning measurements made by the ever-increasing use as reference sources 
of artificial satellites or extra-terrestrial objects such as quasars. The modern 
measurement types include: very long baseline interferometry (VLBI); lunar and 
satellite laser ranging (LLR and SLR respectively); and "ranging" to the Global 
Positioning System (GPS) of satellites. In all these techniques, the observed objects are 
used to approximate a non-rotating reference frame either directly, as in the case of the 
slow moving quasars, or from the dynamical theories of their motion in the case of the 
moon and the artificial satellites. Notably, VLBI and SLR have the capability of providing 
three-dimensional coordinates of relative positions of reference points with an accuracy 
of 5 to 10 cm or better, practically for any baseline length. Repeated observations at 
that level of accuracy can indicate variations in the positions due to crustal dynamic 
processes of the earth on both regional and global scales, and if carried out at frequent 
time intervals, variations in the earth's rotation. 

Fixed VLBI or SLR stations alone can be used to study such motions if they are 
located in the geophysically active regions, and to contribute to the maintenance of a 
reference framework with respect to which the motion of additional points of interest 
can be determined by means of mobile units. In fact, mobile systems are being used 
extensively within the current activities of NASA’s Crustal Dynamics Project (CDP) to 
measure baselines between the fixed sites and locations accessible to the mobile systems, 
thus resulting in much stronger networks and making it possible to directly assess 
deformations associated with tectonic activity in regions within plate boundaries and 
their peripheries. Such observations are extremely valuable and well suited to the 
capabilities of these techniques. Today, however, there is an increasing need for more 
geodynamic support requiring accuracies of the order of 0.1 to 0.01 ppm over distances 
from 100 to 10000 km (Hothem and Williams, 1985) and at more frequent intervals 
both in a spatial and temporal sense. Such measurements could provide the necessary 
information for the adequate characterization of plate velocities and the complex 
patterns of crustal strain so necessary for the full understanding of the tectonic 
processes occuring at the interplate boundaries. VLBI and SLR networks in this context 
are expensive propositions, however, since rapid repeat measurements are difficult to 
carry out. In the near future, it may be possible to monitor geodetic networks with 
spaceborne laser systems that range to corner retroreflectors on the earth's surface 
(Cohen et al., 1986). The feasibility of such systems, however, has not yet been 
demonstrated and the state of current technology suggests that their deployment may be a 
decade or more away. 

This role may well be fulfilled by using alternative techniques such as those 
utilizing the signals from GPS which, already without the full deployment of the system, 
has the same uncertainty level for baselines less than about 1000 km in length as 
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mobile VLBI or SLR. Given the small, relatively inexpensive and highly portable 
receivers, along with the shorter observing requirements and hence, the smaller 
operating costs of GPS campaigns, it appears that GPS may offer a definite advantage for 
this type of application. A similar impact is expected to be realized in the areas of the 
conventional geodetic networks. In their educated opinion, Vanicek et al. (1984) and 
Delikaraoglou and Steeves (1985) predicted that the shift to relative positioning over 
thousands of kilometers instead of a few tens of kilometers or less, will eventually give 
way to a simple reference framework of accurately positioned sites that will act not only 
as foothold points for the relative positioning techniques, but also provide the link to 
specific coordinate systems. In fact, just such a reference framework is currently being 
tested in both the U.S., as part of the validation of JPL’s proposed Fiducial Network 
System concept (e.g. Thornton et al, 1986), and in Canada (Delikaraoglou et al, 1986) 
as part of the development of a GPS-based Active Control System (ACS)-see also Section 
5. 


In the remainder of this section we review briefly the VLBI and SLR techniques 
vis-^-vis GPS mainly for the purpose of gaining a perspective in their respective 
application areas and for placing the role of GPS in context. We shall not dwell in detail 
into the theory of these techniques as it can be found adequately described, for example, 
in IEEE (1985), AGU (1985) and Wells et al. (1986). Rather, we shall attempt to 
outline the relative merits and disadvantages of each and explain how their 
complimentary application can be beneficial. 


1.2 Capabilities and Limitations of Ruailable Techniques 


1.2.1 Uery Long Baseline Interferometry (ULBI) 


The technique of very long baseline interferometry (VLBI) has been used 
extensively since the late sixties to provide, among other applications, measurements of 
the earth orientation, global and regional crustal motion, deep-space spacecraft 
navigation, and three-dimensional vector baselines up to intercontinental distances 
(Clark et al., 1985; Davidson and Trask, 1985; Herring, 1986). Using the principles 
of wave interference and bandwidth synthesis, signals from distant quasars or other 
extragalactic sources received in the antennas of two or more radio telescopes are 
amplified and translated to a lower frequency band under control of a hydrogen maser 
frequency standard; in turn they are digitized, time-tagged and recorded on wide 
bandwidth magnetic tapes. Tape recordings from different sites are then played back and 
cross-correlated to obtain estimates of group delay (i.e. the difference in arrival times 
of the quasar signal wavefronts at the radio telescopes) and phase-delay of the sample 
cross-correlation function and its time rate of change or delay rate. Group delays are 
measured in two widely separated frequency bands in order to derive a new observable 
which is free of ionospheric (or more generally, plasma) delay. For geodetic VLBI these 
frequencies are in the microwave portion of the electromagnetic spectrum; e.g. in the S- 
band (= 2.2 GHz) and X-band (= 8.4 GHz) used in the Mark-Ill system. These 
observables depend among other factors on the direction of the source(s) and the vector 
separation and orientation of the antennas of the observing sites. Hence, it is possible to 
estimate these and other quantities by measuring the delays and their rates of change for 
many different sources and many baselines observing simultaneously through a multi- 
parameter least-squares fitting procedure in which geophysical, astrometric and station 
coordinates are adjustable and/or recoverable parameters. 
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To date, VLBI can provide the highest possible geodetic control needed for the 
precise monitoring of crustal deformations. This is mainly due not only to its high 
accuracy but also to its excellent repeatability. The latter is mainly ensured because of 
the very geometric nature of VLBI and the coordinate system defined by the observed 
radio sources whose positions as viewed from an earth observer define essentially a 
quasi-inertial reference system (which is as close as we can get to an inertial system 
once the small undetected proper motions and the earth’s own inherent acceleration in 
space are neglected). In the analysis of VLBI observations, this system is eventually 
related to a geocentric, earth-fixed coordinate system with the orientation of the earth 
being defined by the pole position, UT1, and by the standard formulas for sidereal time, 
precession, nutation and earth tides, etc. In this context there seems to be a generally 
accepted consensus that VLBI is the only technique that would allow the unambiguous 
establishment of the future conventional terrestrial reference frame through its 
realization by a global network of VLBI stations and the global tectonic motions of those 
positions. Satellite techniques like SLR and GPS do not have direct access to such a 
quasi-inertial reference frame. Besides the lack of direct connections of the coordinate 
system actually employed and the inertial frame, these techniques are sensitive by 
contrast to gravitational forces hence, insofar as this knowledge is incomplete, errors in 
computing the force field in which the satellites move could and do lead to incorrect 
satellite positions which, in turn, affect the geodetic coordinates of the ground stations. 
Furthermore, such techniques basically observe range or range differences and contain 
no information about the rotational or translational degrees of freedom of the satellite 
networks, thus resulting in additional estimability problems (Grafarend and Mtiller, 
1985; Delikaraoglou,1985) which can be overcome only after carefully chosen 
geometric or dynamic constraints are imposed. 

The use of VLBI is not limited to large fixed radio telescopes. Since the early 
1980’s, mobile systems have been developed by JPL for the NASA Crustal Dynamics 
Project and have been operated in conjunction, for example, with several fixed large 
antennas in the Western United States and Canada to establish relative motions near the 
tectonic plate boundaries in California and Alaska. The small (from 3.0 m to 9.0 m 
diameter) antennas used with such systems are a major advantage for geodynamic 
applications, particularly in allowing the establishment of control within the global 
networks of the fixed sites and more frequent measurements (both in space and time) to 
adequately determine the deformations within the plates and their peripheries. Like 
their fixed counterparts, mobile VLBI systems have essentially all-weather capability 
but involve considerably more complex operations and poorer system consistency, 
mainly because the smaller size antennas yield less sensitive interferometers than 
normally achievable with the larger fixed antennas. A direct consequence of this is the 
reduction of quasar signal-to-noise ratio which, in turn, imposes severe limitations in 
the observing schedules (often restricting observations to the stronger sources) and 
tends to distort the experimental geometry and observing strategies. This decrease in 
sensitivity can of course in principle be compensated for by higher-gain antennas, 
lower-noise receivers and multi-observing sessions, and indirectly by differencing the 
vectors between individual mobile units and a common fixed site, but not without the 
expense of ail the attendant complexities and increased costs of operations. Naturally, 
there are other error sources inherent in VLBI which impart temporal variations on the 
baseline components. Many of these sources of error are in fact shared by systems like 
GPS as will be discussed in Section 2. Given the complexity of the VLBI operations, costs 
in maintaining VLBI networks and correlator facilities are obviously high, and rapid 
repeat measurements are difficult. This is, however, a small price to pay in light of the 
usefulness of the technique in providing a stable reference system for other geodetic and 
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geodynamic work including the determination of accurate orbits for systems like GPS, as 
will be seen in some detail in Section 5. 


1.2.2 Satellite Laser Ranging (SLR) 


Satellite Laser Ranging (SLR) systems have been in use since the mid-1960’s. 
By tracking earth-orbiting satellites from a global network of laser tracking stations, 
information on the forces acting on the satellite can be derived from the analysis of 
measurements of the round-trip travel time (and hence the path length) of a photon 
pulse from a laser site to the corner-cube reflectors mounted on the surface of a 
satellite. Overall laser tracking system capabilities have increased in accuracy by a 
factor of 3 every 5 years over the last 20 years as noted by Pearlman (1982). To date, 
the distance from a tracking site to the satellite can be determined to a precision of a few 
centimeters using current technology and techniques. Most SLR systems now use short 
pulse length (100-200 ps) neodymium YAG lasers with receiving sensitivities ranging 
from a single photon to more than 100 photons and best single-shot range accuracy at 
the level of 2 cm. 

Currently there are some 25 fixed SLR systems worldwide. In addition, NASA 
has developed and deployed eight trailer-based mobile Laser systems (MOBLAS) and four 
highly compact transportable-class (TLRS) systems with a further two transportable 
sytems (MTLRS-1 and -2) having been built in Europe. The TLRS systems were 
developed mainly in response to a need of obtaining SLR data at remote sites and are 
rapidly being relocated in support of local and regional tectonic studies. Much like the 
mobile VLBI systems, the compact size of the TLRS systems was achieved at the expense 
of greatly reduced operational signal levels (about two orders of magnitude lower) as a 
result of the smaller telescope apertures (typically 25 cm vs. 75 cm) and lower laser 
output energies (5mJ vs. 100 mJ) relative to the larger MOBLAS systems (Degnan, 
1985). 


Of the retroreflector-equipped satellites the most fruitful in terms of geodetic 
results, have been the dedicated, fully reflecting laser satellites LAGEOS and STARLETTE, 
and only recently the Japanese Experimental Geodetic Payload satellite, AJISAI. Both the 
geodetic and geodynamic inferences drawn from SLR data depend mainly on dynamic 
techniques which require accurate descriptions of the earth's gravity field and the 
precise definition of the coordinate frame of the tracking network locations, including 
the description of polar motion and specification or determination of fundamental 
constants (e.g., GM, speed of light, etc). In fact, for current laser tracking systems, 
errors in the gravity field models are one of the dominant errors in the recovery of site 
locations and certain orbital and other geophysical parameters. Several other error 
sources are also present that limit the accuracy of the present systems. These include 
uncertainties in the relative timing errors, atmospheric refraction, epoch time errors, 
system delay calibration errors, etc. There are also other factors that are detrimental 
to SLR activities. Although susceptibility to adverse weather is a major drawback, 
system availability and scheduling are also often cited as limiting the efficiency factor at 
about 30% for the current systems. As a result, operational costs are particularly high 
for SLR. Typically, 5 to 30 days of occupation for fixed sites (and up to 60 days for 
mobile systems) are required to determine the length of baseline vectors up to 
intercontinental distances to a precision level of 3-5 cm. 
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In spite of these disadvantages, SLR systems have provided very precise 
measurements for establishing the orientation of the earth within a quasi-inertial 
coordinate frame (Smith et al., 1987). This is ensured in part, by the methodology used 
which interrelates the laser observations through a theory of motion describing the 
satellite’s orbital behaviour in time and subsequently determines the positions of the 
tracking stations in an appropriate terrestrial coordinate system. The latter, being on 
the earth’s surface, rotates with it while the satellite orbit plane is, to a first 
approximation, independent of the rotation of the earth. Hence, the satellite motion 
(described in the inertial frame), and the stations (located in the earth-fixed coordinate 
frame) are linked through the earth’s angular velocity vector which is oriented by the 
polar motion and UT1 solutions. In the case of SLR, particularly for the LAGEOS 
satellite, the data are sufficiently abundant to actually permit a complete solution of 
station positions and earth orientation parameters simultaneously with the satellite 
orbits-typically at a level of precision of 2 mas for the x,y axis of the earth’s rotation 
axis and 0.2 msec for UT1 for data in 5-day intervals (Tapley et al., 1985). 


1.2.3 Global Positioning System (GPS) 


The Global Positioning System (GPS) of satellites is currently being developed by 
the U.S. Department of Defense (DoD) as the next navigation system of the future. A 
detailed description of GPS may be found in ION (1980, 1984), NOAA (1985), ARL 
(1986) and Wells et al (1986). In summary, the system is designed to provide all- 
weather, real-time precise position and time information with a higher degree of 
accuracy than ever before available. When fully operational (in the early 1990’s as 
now postulated following the delays of the launch schedule caused by the Space Shuttle 
Challenger accident), the satellite segment of the system will consist of 18 operational 

satellites and three active spares placed in six orbital planes inclined at 55° to the 
equator (the present prototype satellites are at inclinations of 63° instead) at a height of 
about 26,000 km, having orbital periods approximately 12-hours, thus ensuring that 
any place on the earth will be able to "see" a minimum of four and up to eight satellites 
at any instant in time. 

The GPS satellites carry an ensemble of atomic clocks onboard which are 
carefully monitored by the U.S. Naval Office (USNO) and their time scales are 
maintained trackable with respect to UTC-USNO to within ± 100 nsec, which ensures 
high satellite oscillator stability and precise timing information. Each satellite 
transmits at two L-band frequencies, Li (1575.42 MHz) and L 2 (1227.60 MHz) which 

allows for the elimination of the first-order ionospheric effects. Three different spread 
spectrum codes are modulated onto the L-band signals, including a data (D-code) 
modulation consisting of a satellite broadcast ephemeris message (at the rate of 50 
bits/sec) and two pseudo-random noise (PRN) codes known as the coarse/acquisition 
(C/A) and precise (P) codes. The latter are modulated onto the L-bands at frequencies of 
1.023 MHz and 10.23 MHz respectively, with the respective codes repeating every 1 
msec for the C/A code and approximately every 38 weeks for the P code (although in 
practice the P code for each satellite is reset every week). It is the current DoD policy 
that the future availability of these codes, when the operational system is in place, would 
eventually rely on the so-called Precise Positioning Service (PPS) to be based on the P 
code on two frequencies and the Standard Positioning Service (SPS) to be based on the 
C/A code only on the L 1 frequency. In order for the SPS to be freely available, it will 
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necessarily be poorer in accuracy, while the PPS will be more precise but will be 
encrypted and thus limited in access to civil users (Baker, 1986). Furthermore, it is 
entirely conceivable that DoD may severely degrade the relative accuracy of the SPS in 
the interest of national defense. 

The present (1987) status of the system is that 11 (prototype) satellites are 
currently in orbit, with 7 satellites being fully operational allowing 4 to 5 hours per 
day of four-satellite coverage in the United States and Canada, Europe, a large portion of 
South America, Australia, the Near East and the polar regions. Even though the system is 
still technically in its testing phase, it has already demonstrated its capability for 
determining baselines with an accuracy of 0.1 ppm over distances of a few hundred 
kilometers and, most importantly, with much shorter required observing periods for 
such accuracies than with other techniques. This, together with the all-weather 
capability and the true portability of the currently available GPS equipment combined 
with the relatively low cost of operations, clearly illustrate the apparent advantages of 
GPS. This should not, however, mitigate the few disadvantages inherent in the system as 
well, particularly the dependence of the user on the policy of the DoD. Possible 
arbitrary, sudden system shutdowns and/or failures and lack of precise orbit 
ephemerides are extreme examples of this policy which are obviously nonexistent in 
systems like VLBI and SLR. The possibilities of selective availability and degraded 
accuracy should have some effect on static relative positioning (especially on 
establishing appropriate modelling strategies) with the major factors affecting the 
accuracy of relative coordinate determinations being: 

- the accuracy of the orbit representation, and 

- our capability to adequately model atmospheric refraction errors, particularly 
those of the wet troposphere and especially over long baselines, 

which would be very critical in the context of the high accuracies needed for geophysical 
and geodynamic applications. In spite of these drawbacks, however, GPS remains an 
extremely attractive measuring technique likely to fulfill the future role of mobile VLBI 
and SLR to provide high-precision densification control in the local and regional areas of 
the VLBI and SLR networks. This premise has led the U.S. National Geodetic Survey 
(NGS), NASA, JPL, and others to explore this feasibility through several campaigns that 
have been initiated to use GPS for the study of crustal motion (e.g. Strange, 1985; 
Thornton et al., 1986). In Section 5 we shall report on some of the most recent results 
of GPS relative positioning in terms of the levels of accuracy achieved by employing the 
methodology described in this report for minimizing the various error biases, 
particularly those of the satellite orbits. 
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SECTION 2 


THE RCCURRCV OF GPS RELRTIUE POSITIONING 


Since GPS has been designed primarily as a military navigation system to provide 
all-weather, real-time navigation, its principal mode of operation is real-time point 
positioning by ranging using the codes of either the available PPS or SPS signals. This is 
the way GPS is designed to be used for navigation, although even in this mode it could 
provide geodetic accuracies if used on the same station over an extended period of time 
(Wells et al., 1986). The type of observable, however, which can be used to perform 
precise relative positioning independent of the coded signals is that of the carrier beat 
phase. Carrier beat phase measurements with GPS are similar to the continuously 
integrated Doppler count measurements obtained from the Doppler TRANSIT satellite 
system. That is, a carrier phase measurement is essentially the phase of the signal 
which remains when the incoming Doppler-shifted satellite carrier signal is beat with 
the nominally-constant reference frequency generated in the receiver. Clearly, if the 
transmitter remains at a fixed distance from the receiver , the so-measured phase of the 
incoming signal remains constant. Viewed from this perspective, it is evident that the 
phase of the signal changes with the change of the range between the satellite and 
receiver: the phase changes by 1 cycle wavelength of the carrier frequency (19 cm on 
the Lf and 25 cm on the Lg GPS frequencies respectively). This type of measurement 

leads potentially to the most precise information about the satellite-to-receiver ranges 
one can obtain from GPS. However, the problem of utilizing this potential is one of 
ambiguity: it is quite difficult to accurately locate the cycle of the carrier whose phase 
is being measured. With currently available instrumentation, phase changes measured 
to an accuracy of a fraction of a cycle are possible. Hence, carrier beat phase 
measurements are sensitive to sub-decimeter range changes. However, the success of 
achieving this level of positioning accuracy from this type of data hinges on the 
capability to resolve this cycle indeterminacy (or ambiguity) including the recovery of 
phase cycle slips often encountered during tracking. 

Relative or differential GPS positioning requires that two or more receivers are 
observing simultaneously the same satellites. Given the altitude of the satellites in the 
system, this implies a limit of about 4000-5000 kilometers on the length of baselines 
that can be measured and also that the data on both stations must be correlated for a 
solution. One can envision various levels of such differential operations. At the lowest 
level, it is necessary only that the operating receivers track GPS at the same time, 
although they need not track the same satellites. In this case, the effect of the GPS timing 
and propagation errors would be correlated. One level higher, the receivers would track 
the same satellites, although not necessarily in the same sequence. In this case, the 
effects of ephemeris errors and specific satellite timing biases would be correlated. At 
the highest level, the receivers would be forced to track the same signal wavefront from 
the same satellite simultaneously, in which case the correlation between the effects of 
satellite and propagation error sources would be at a maximum. This latter mode of 
operation, which is very similar to the conventional VLBI, is most often used for geodetic 
and geodynamic applications, since it leads to one or two orders of magnitude 
improvement in the results. 

Differential observations can be processed on an individual baseline-by-baseline 
basis or all together as a network. The data processing can furthermore be done either 
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with the direct "raw" phase observable (i.e., in an undifferenced mode) or with some 
linear combination of the measured phases from the different stations. The linear 
combinations are usually obtained from the raw phases which are differenced further 
between-stations, between-satellites, and between-epochs and in combinations thereof 
that remove, or greatly reduce specific types of errors. A single difference is the 
difference in ranges (or phases) from two receivers observing simultaneously the same 
satellite. These observables cancel any satellite-specific biases such as satellite clock 
offsets and to a great extent (depending on baseline lengths which are short compared to 
the altitude of the satellites) satellite ephemeris and atmospheric refraction errors. 
Double differences are the differences between two single differences of two satellites. 
These observables tend to further cancel any receiver-specific biases such as receiver 
clock offsets (which are common contributors to both single differences), but they are 
mathematically correlated with others in the same stations or to the same satellites. One 
differencing step further, triple differences are obtained by differencing between double 
differences over two consecutive time epochs. This combination also cancels any 
residual satellite and receiver clock errors and any receiver-satellite pair-specific 
biases such as cycle ambiguities. Like the double-difference observables, triple 
differences are also mathematically correlated. 

Theoretically, subject to including appropriate covariance matrices reflecting 
these correlations, the data reduction results with any of these observables should be 
identical. Conversely, neglecting these correlations could lead to significantly different 
results (Ashkenazi and Yau, 1986). Undifferenced carrier phases may also be processed 
directly by including extra nuisance parameters in the model to account for the various 
biases (Goad, 1985; Lindlohr and Wells, 1985). Furthermore, it has been shown that 
the undifferenced and differenced techniques are equivalent, if the same biases that are 
cancelled in the differenced observables are included as nuisance parameters in the 
undifferenced model (ibid; Bock et al., 1986; Schaffrin and Grafarend, 1986). In 
practice, however, these methods will give different results, particularly over the 
longer baselines (> 200-300 km) which are of interest to geodynamics, mainly because 
significant atmospheric variations reduce the correlations among the observations, thus 
offsetting the advantage of using such observations to reduce the existing errors in the 
first place. In fact, this is the main reason why the use of double differences weighted 
according to baseline length is argued (e.g. by Bock et al., 1986) to be the preferred 
method of processing since the residual atmospheric effects in these observables are 
approximately proportional to baseline length. 

As the discussion has alluded to already, there are two major types of 
requirements in GPS operations that differential techniques are capable of supporting. 
Firstly, differential GPS must be capable of improving the performance of the system 
beyond the limits of its non-differential expected accuracy envelope. This is brought 
about merely by utilizing the differential observables just mentioned. Secondly, 
differential GPS must merely be capable of detecting when GPS is outside its expected 
accuracy range, whether this is due to intentional or unintentional degradation, signal 
propagation, timing, orbital or geometric causes, all of which gain particular 
significance in precise operations for geodynamics. This requires particular attention 
in understanding the "cause and effect" relationship between the various biases and 
errors that are inherent in the GPS measurements. 

In the remainder of this section, we shall attempt to identify and describe, mostly 
in non-mathematical terms, these error sources with respect to the accuracy of their 
modelling and to describe current efforts to deal with their effect in general, and in the 
present study in particular. 
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2.1 Biases and Errors 


The accuracy of positions obtained by GPS is dependent on two general influences: 
the errors affecting the measurements themselves and the geometric strength of the 
relative satellite-station configuration being observed. The geometrical aspects relating 
to the latter influence have been discussed, e.g. by Vanicek et al. (1984). Hence, nothing 
more will be said about them here since they are of no direct consequence to the present 
study. Later we shall only discuss the relation between orbital errors and relative 
position accuracies as it can be established from a purely geometric point of view. 

The measurement errors generally fall into three categories: those associated 
with the satellites, with the signal propagation, and with the receivers. Errors from 
each of these sources have complex spectral properties and there are correlations 
between some of them, both in the temporal and spatial sense. These will be discussed 
immediately below. 


2.1.1 Satellite Orbital Biases 


Satellite-associated biases consist of errors in the orbital ephemerides (i.e., the 
satellite is not where the GPS ephemeris used in the data reduction tells us it is) and 
errors in the satellite clocks (i.e., the satellite time scale is not perfectly synchronized 
with the various other time scales usually involved in the data reduction-see also 
Section 4). These errors are uncorrelated between satellites; they affect both code and 
carrier beat phase measurements equally; and they depend on the number and the 
location of the tracking stations providing data for orbit determination and monitoring of 
the satellite clocks (i.e., the GPS control stations), the force models used and the 
satellite configuration geometry itself (Swift, 1985; Fliegel et al., 1985). Satellite 
ephemeris errors are the most difficult to deal with. Clocks would be improved, or their 
undesirable effects would be eliminated or greatly removed by the differencing of 
simultaneous observables from two or more receivers and satellites. Ephemeris errors 
on the other hand, would require better understanding and estimation of the forces acting 
on the satellites; a process which is hampered however, by insufficient knowledge since 
these forces cannot measured directly from the ground tracking sites. The main 
approach to solving the problems created by the imperfect knowledge of the physical 
phenomena responsible for these effects has been usually directed towards developing 
more “exact” models to begin with. This is the standard approach followed by many to 
date. In this study, while following the same example in making the same choice, we 
also have opted for the alternative approach (cf. Section 3): to use parametric models of 
various complexities for the ephemeris parameters themselves (i.e., parameters that 
can be adjusted as part of the orbit estimation process itself). 

The broadcast ephemeris is known to be in error, typically by 20 to 80 m, which 
clearly makes it unsuitable for any precise applications, particularly over long 
baselines. It is also now known that deficiencies exist in the solar radiation force model 
presently employed in the estimation and prediction of the orbits of the prototype 
satellites by the Naval Surface Weapons Center (NSWC), where the so-called precise 
ephemeris is also being computed, and that a new improved model is to be required for 
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the operational satellites (Fliegel et al., 1985). Although the accuracy of the NSWC 
precise ephemeris is not exactly known, in an earlier investigation, Swift’s (1985) 
evaluation indicated that it may be in error by as much as 6 m, 43 m and 25 m in the 
radial, along-track, and across-track directions respectively. As soon as the 
operational configuration of the satellites is in place, it is expected that the operational 
ephemeris accuracy will be of the order of 1 .5 m in the radial direction (Van 
Dierendonck et al., 1980), although this is still thought to be quite optimistic given the 
current experiences. Rather, actual errors of about 5 m are thought to be more likely 
(Wells et al., 1986). 

The relation between orbit errors and resulting relative positioning accuracies 
can be established from a purely geometrical point of view. The exact linear relation 
between the basic differential range observable ApJ and the baseline vector B between 
two sites can be expressed as shown in Fig. 2.1 (cf. also Vanicek et al., 1984) by 


ui B = - ui eJ 1 ApJ 


- - [1 - 14 CO 4 ] Apj 

where 


ui = ^(e^ + ei 2 ) 


( 2 . 1 ) 


( 2 . 2 ) 


is the average of the unit vectors ei-j and ei 2 in the direction of the satellite. 

Assuming that an error in the satellite position introduces an error into the 
baseline vector in such a way that the observable is unaffected, one can obtain a first- 
order approximation of the error equation that relates these errors by replacing ei and 
B by eJ+dei and B+dB. This results in the error equation 


(dei 1 + dei 2 ) B + (ei 1 + eJ 2 ) dB - 0 


(2.3) 


which defines a constraint for the effect of the orbital errors on a baseline determination 
and makes a unique solution for dB possible if the condition that the magnitude of dB be a 
minimum is imposed along with (2.3). The form of this error equation indicates that 
the errors introduced into the baseline vector depend on the magnitude and the 
orientation of the orbit error vector relative to the baseline vector B . Buffett (1 985) 
has shown that the solution of the error equation that maximizes the effect of the orbit 

errors is obtained when dB is parallel to ei^ and eJ 2 which corresponds to a more 
convenient rule of thumb that has also been adopted by several investigators (e.g. 
Vanicek et al., 1985), i.e. 


(dB/B) = (dp/p) 


(2.4) 
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Fig. 2.1 - Differential range geometry. 


where dp is the orbit error expressed in terms of a topocentric position error and p is 
the average topocentric distance (=20,000 km) to the satellite. The implication of the 
rule of thumb for relative positioning is that although for short baselines (< 200 km) 
sensitivity to the orbit errors is reduced by a factor of 100 more, for centimeter 
geodesy over distances greater than a few hundred kilometers, orbit errors still limit 
system performance. This implies that for regional and global geodynamic applications 
which require 0.01 ppm accuracies (i.e., 1 cm over 100 to 1000 km), the satellite 
orbit errors should not exceed 2 m (in equivalent range error) for centimeter 
accuracies over 100 km and 20 cm for the same accuracy over 1000 km. Consequently, 
the need to improve the GPS orbits has been recognized early on to be a necessary part of 
improving the current GPS capabilities. 


2.1.2 Propagation Biases 


Propagation errors consist of refraction errors due to the influence of the 
ionosphere and the troposphere, multipath, imaging and phase-center variations due to 
the mutual interference of two or more signals emitted from the same source and 
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reflected along different paths by surrounding reflective surfaces in the vicinity of the 
GPS antennas. 


2.I.2.1 Ionospheric delay errors 


GPS code measurements are subject to ionospheric group delay, and carrier beat 
phase measurements to ionospheric phase delay. At the GPS frequencies, the plasma of 
the ionosphere (50-1000 km above the earth surface) ionized by ultraviolet radiation 
from the sun can result in range delays varying from 50 m up to 150 m (at maximum 
sunspot activity, midday, satellite near horizon) to less than a meter (sunspot 
minimum, night, satellite at zenith). Since ionospheric refraction is frequency- 
dependent, measurements on the Li and L 2 frequencies can be compared to estimate the 

first-order effect represented by an expression of the form (Spilker, 1980) 


a/f 2 + b/f 3 + c/f 4 + ... 


(2.5) 


where f is the carrier frequency and a,b,c,... are parameters which are functions of the 
total electron content along the signal path. For the carrier beat phase measurements, 
for instance, this correction to the Li phase observations takes the form 


d<(>i = f 2 2/ ( h 2 ' U 2 ) (2.6) 


where <J>i , <J> 2 are the phase observations in the Li and L 2 frequencies. This first-order 
correction removes about 99% of the range delay effect, leaving some 8 cm at the most 
as a result of neglecting the higher order effects b/f 3 and c/f 4 (ibid). However, there is 
an undesirable by-product of this correction: noise on each of the Li and L 2 phase 

observables will result in increased measurement noise in the corrected observation. 
Assuming for the L 2 a noise amplification of 1.35 times that for the Li observable, a 

simple propagation of errors using (2.6) reveals that the noise level of the dual- 
frequency estimate is 3.3 times that of the uncorrected Li observable which illustrates 

that there is a point at which making the dual-frequency correction becomes 
counterproductive. In differential positioning the corresponding effects at the ends of a 
baseline are obviously of concern. Clearly, for short baselines (<50 km) where there is 
very high correlation of the effect, ionospheric errors can be greatly reduced simply by 
single differencing or even be ignored altogether. This has been verified experimentally 
by several investigators (e.g. Kleusberg, 1986). Of course for double-difference 
observables, there is a further amplification of the noise level due to the linear 
combination for both the dual-frequency correction and the double-differencing 
operation. Thus, the use of dual-frequency combinations is generally not advisable for 
such short baselines. For longer baselines of 100 to 1000 km, there is a different 
situation since signals between satellites spaced far apart (for optimal geometry) will 
generally pass through different and uncorrelated parts of the ionosphere: thus one is 
obliged to remove the ionospheric effect using the dual-frequency correction. In this 
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situation the horizontal gradients of the ionosphere will be the determining factor of the 
residual effect. Usually, the largest gradients will be in the north-south direction at 
mid-afternoon and in the east-west direction at sunrise, and will reach values between 
about 4 m/1000 km during the maximum solar activity cycle and 1 m/1000 km during 
solar minimum (Clynch and Coco, 1986). 

A simple way to account for these residual effects is to consider that the total 
ionospheric effect may be expressed as a function of the total electron content N e of the 

ionospheric layer as 



dB 


Fig. 2.2 - Influence of ionospheric refraction on GPS observations. 
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8 P ion = (k N e / f L 2 ) « (k N v / f L 2 sin E') + (k 8N V / f L 2 sin E') (2.7a) 


- IER + (IRE V / sin E’) - IER + IRE S 


(2.7b) 


where k is a known constant (k=40.5284), fL is the carrier frequency in hertz, N v is 
the estimate of the vertical electron content and 5N V denotes an error in N v . IER is the 
first-order effect accounted for by the two-frequency correction (2.6), while IRE S 
expresses the residual ionospheric error in terms of the vertical electron density error 
IRE V . The simple mapping function 


IRE S - IRE V / sin E' 


(2.8a) 


is a result of assuming a thin ionospheric layer at a mean ionospheric height (usually 
set at H * 400 km; i.e. 50 km above the layer of maximum density in the mid-latitude 
regions) onto which the entire electron content is assumed to be condensed; E' is the 
elevation angle of the slant range at the crossing point of the ionospheric layer. In more 
rigorous terms, the ionospheric crossing point should be described by two parameters: 
the station latitude and the local solar time (given that the ionosphere is strongly 
dependent on the hour angle of the sun; i.e. day vs. night observations). In simple terms, 
however, it can be determined as a function of the satellite elevation and the geocentric 
station vector R as 


cos E' = [R / (R+H)] cos E 


(2.9a) 


The residual ionospheric error IRE S can then be computed as 


IRE S = RREN X (IER) 2 


(2.8b) 


where RREN is a normalizing factor approximated as 


RREN = a + b e- EV20 - 7 


( 2 . 10 ) 


where a, b are coefficients to be determined in the course of the adjustment of the data 
along with other nuisance parameters. 

Unmodelled ionospheric effects in GPS carrier beat phase observations lead to an 
increased noise level by a factor of 2-3 in the adjustment residuals (cf. Kleusberg, 
1986) and to a scale effect on the so-determined baselines (cf. Beutler et al., 1987). 
As a consequence, a baseline B' derived from single-frequency GPS observations could 
be shorter with respect to the true baseline B as shown in Fig. 2.2 by an amount (dB) 
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proportional to the difference of the ionospheric error (dpj on ) 12 = 5p 2 k - Sp^ at the 
two ends of the baseline; i.e. 


dB cos dA = (dpj on ) 12 / cos E 


( 2 . 11 ) 


where dA and E are the difference in the corresponding azimuths and the average 
elevation angle to the satellite respectively. The combined effect of all observations can 
be simply determined by the geometrical configuration and the exhibiting balance of all 
azimuths and elevations over all satellites observed within a given session; i.e. as a rule 
of thumb, 


AB/B = £ l(dpj on ) 1 2 / B cos dA cos E] { (2.12) 

i 

expressed in mm/km. This amounts to about 0.25 ppm per 1 m of unmodelled (or 
improperly modelled) ionospheric error. For numerical computations, Georgiadou and 

Kleusberg (1988) have shown that (dp j on )-i 2 can be conveniently expressed as 


(dpj on )i 2 = dV B [R/(R+H) 2 ] sinE cosE (cosz')‘ 3 x 

x {1 + R/r sinE + (R/r sinE) 2 } (2.13) 

where r is the geocentric distance of the satellite, z' is the zenith distance corresponding 
to the given elevation angle E and computed from (2.9a) as 

sin z' = [R/(R+H)] cos E ; (2.9b) 

the vertical ionospheric delay dV is a function of the electron content 

dV - 40.3 N v / f L 2 (2.14) 


which itself is the most influencing factor given that its variability is the strongest both 
spatially and temporally as well as regionally. Layer height plays an important role, 

particularly in elevations of about 15°. For a vertical delay of 1 m and a baseline of 1 
km, (dpj on )^ 2 as computed from (2.13) may vary in the fashion illustrated in Fig. 2.3 

showing the dependence of the effect on the assumed height of the ionospheric layer at 
three altitudes: H=300 km (low), H=400 km (normally assumed) and H=500 km. 
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Fig. 2.3 - Error in differential ionospheric path delay 
per 1 m of unmodelled (or improperly modelled) vertical ionospheric delay. 


At sites where single-frequency GPS receivers might be operating, the vertical 
ionospheric delay, dV, along the paths of the observed satellites can be determined from 
the carrier phase observations of one or more dual-frequency receivers located in the 
vicinity of the GPS survey. This is a practical application where the usefulness of an 
Active Control System of GPS sites can be indeed very significant. In fact, some 
experimental analyses of data collected in Canada over short baselines (< 30 km) have 
already been conducted (ibid) to test this premise. Although positive indications were 
drawn from the results of these tests, further tests are necessary over long baselines 
before any definite conclusions are drawn on this issue. 


2. 1.2.2 Tropospheric delay errors 


GPS signals, like VLBI, are experiencing a delay due to the neutral atmosphere 
(which includes the troposphere and other regions up to some 80 km altitude) which is 
essentially independent of the frequency of the carrier. The tropospheric contribution of 
the slant range to a satellite is usually expressed as 


(Sp)trop = (^P)dry ^(^)dry + (^P)wet ^(^)wet 


1 6 


(2.15) 




where (5p) dry and (8p) wet are the contributions of the dry and wet troposphere at 
local zenith and R(E)<j ry and R(E) we t are corresponding elevation angle mapping 
functions. The modelling of the effect of the dry troposphere is easier and fairly 
accurate to handle as demonstrated by SLR and VLBI data analyses. On the other hand, the 
wet portion of the troposphere is much more difficult to model mainly due to 
uncertainties in determining the representative humidity along the signal path which is 
not necessarily well correlated with surface conditions. This error can be reduced by 
using water vapour radiometers (WVR) which, unfortunately, are expensive and may 
not operate under all weather conditions or offer direct external calibrations of the low 
elevation angle observations necessary in continental scale GPS experiments due to the 

low elevation limit of 20° for current instruments. 

The best models available to date can typically reduce the combined effect by 92- 
95% or leaving up to 2-5 cm on the average, depending on the amount of atmospheric 
information available. In more extreme cases, e.g., at typical humid coastal sites, Ware 
et al. (1987) have found, for instance, 12-cm rms errors in the vertical component of 
GPS baselines when surface measurements were used in a tropospheric model similar to 



Baseline length (km) 


Fig. 2.4 - Effect of wet tropospheric refraction modelling 
errors on GPS results translated from VLBI experiences based 
on observation periods of 24 hours (after Kouba, 1987). 
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(2.15) to estimate path delay. In most GPS data reductions, zenith tropospheric delay 
parameters are statistically estimated. This technique essentially allows accounting for 
a spatially and temporally average troposphere for each station and for the time over 
which the zenith scale parameters are estimated. The wet tropospheric variations 
around these averages are then the dominant tropospheric errors which, if unaccounted 
for, in turn map into errors in the baseline length and height differences between 
observing sites. For relative positioning GPS operations over baselines less than 200- 
300 km other errors such as orbital and ionospheric errors have often obscured the 
influence of the wet tropospheric mis-modelling on the determined baselines. 
Translating the VLBI experience into the GPS situation, Kouba (1987) has suggested a 
simplified description of the effect of the wet troposphere on the determined baseline 
lengths and height differences between GPS sites using a model of the form 


a 2 j = s 2 [l-exp(-B 2 /d 2 )] / B 2 


(2.16) 


where (jj is in mm/km, the index i stands for baseline length or height difference, and B 
is the baseline length in km. Equation (2.16) assumes that the effect is governed by a 

second-order Gauss-Markov process described by a discrete variance s 2 (s = 40 mm for 
length, s - 80 mm for height differences which correspond to long baseline uncorrelated 
tropospheric conditions based on VLBI experience) and a correlation distance d usually 
set at 30 km which (although somewhat arbitrary) reflects the common portion of the 
wet atmosphere seen from the two ends of the baseline. The underlying property of this 
process is that the baseline correlation diminishes with increasing distance and 
approaches unity as the station separation approaches that of a zero baseline. Numerical 
values from the effect computed from Eq. (2.16) are plotted in Figure 2.4 showing that 
short baselines are subject to a systematic error of about 1 ppm for baseline length and 
nearly 3 ppm for height differences. For distances over 50 km the effect decreases 
proportionally to the inverse of baseline length, vanishing at about 500 km which 
corresponds to the VLBI results. At distances greater than 500 km the wet tropospheric 
scale and height effects are assumed to be random under the assumptions of the 
underlying Gauss-Markov process upon which these values are based. Treuhaft and Lanyi 
(1985) have shown that if a more detailed model of atmospheric structure is used, 
water-vapour fluctuation-induced errors are caused mainly by temporal fluctuations 
due to spatial patterns which are moved over a site by the wind and hence, vary with the 
observing schedule. In Treuhaft’s and Lanyi’s educated opinion such errors can be 
reduced at their optimal level if the length of time for the zenith solutions is less than 
approximately 4 hours, which coincidentally is typical for the GPS observing sessions. 
To date, to the author’s best knowledge, not enough experiments have been conducted to 
assess the character of the noise contribution of the WVR instruments vis-^-vis that of 
the GPS data and to resolve the issue of the usefulness of WVR’s for individual GPS 
experiments. In that light, the practical benefits of using a calculated tropospheric 
covariance according to (2.16) or similar models for the parameter estimation scheme 
should be studied further. 



2.2 Multipath and Imag in g 


Multipath effects result from the mutual interference of two or more signals 
emitted from a satellite that travel along different paths due to surrounding reflective 
surfaces (Fig. 2.5), thus arriving out of phase in the antenna causing interfering signals 
which may not be decodable or understandable by the receiver. Antenna imaging is a 
secondary phenomenon, similar to multipath, caused by nearby conducting objects 
coupling electrically with the antenna, thus essentially defining an "image" of it with the 



Reflective Surface 

Fig. 2.5 - Geometry of multipath effect on GPS signals. 


resulting amplitude and phase characteristics (particularly of the phase center of the 
antenna) being significantly different from those of the isolated antenna. 
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Based among other factors on the geometry of the satellites and the antenna 
characteristics, multipath effects on carrier phase observations exhibit themselves as a 
very irregular function of time. The phase difference between two "out-of-step" 
arriving signals is given directly as (Tranquilla, 1987) 


O = 2rc AL /X 


(2.17) 


in rad/m, where AL is the path length difference. Assuming that the antenna 
characteristics remain unchanged, the interference patterns of direct and reflected 
multipath waves will exhibit themselves in the data time series as we move through the 
observation plane. For instance, for a series of carrier phase GPS observations on a 

single satellite (cf. Fig. 2.5) corresponding to an observation angle sweep of 5° (or a 
time interval of 5-10 minutes) say, at u-|=45 0 and U2=50°, the multipath-induced 

phase change is computed from Eq. (2.17) approximately as 0.064 x 2 k AL/X rad/m 
which in turn, after setting it equal to one period implies AL/X=15.5 or approximately 
AL=3 m at the Li frequency; clearly a very significant effect. In practical terms, this 
corresponds to a single reflective surface located at a distance of 3 m. A reflective 
surface at a closer distance would produce low-frequency variations, whereas a surface 
at a greater distance would produce high-frequency variations. This has been confirmed 
by experimental observations (e.g Bletzacker, 1985; Evans, 1986; Tranquilla et al., 
1987) which have shown that multipath effects can induce periodic effects ranging from 
sub-minute to 2-3 min periods (caused, for example, by diffused forward scattering 
multipath such as from overhead wires at low elevation angles); 5-10 min periods 
varying with satellite elevation (caused, for example, by specular reflections from 
reflective surfaces in the vicinity of the antenna); or up to 25-60 min periods (usually 
associated, for example, with reflections from fresh or saltwater surfaces). 

Although it has been suggested that the only means to minimize these undesirable 
effects is through antenna beam shaping (e.g., the use of polarized antennas, large ground 
planes and sheets of absorbing material), it may be possible to reduce the effect of 
multipath either by averaging over long enough time periods, so that the relative phase 
of the direct and reflected signals changes by at least one full cycle or alternatively by 
accommodating spectral analysis techniques during preprocessing to detect and suppress 
any periodic variations due to multipath found in the data series. Such an approach has 
been tested in the software developed in the course of this study and seems to work well 
when dealing with this problem at the stage of cycle-detection in the preprocessing of 
carrier phase data. More will be said about this below. 


2.3 C ycle Slips 


As pointed out in the beginning of this section, carrier beat phase measurements 
lead potentially to the most precise information about the slant ranges to a satellite one 
can obtain from GPS. The problem of utilizing this potential however is twofold; 
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(a) Associated with the observations is an ambiguity in the number of whole 
cycles between the satellite and the receiver at the time of initial signal acquisition (see 
also Section 4.3 and Fig. 4.2). Although with presently available instrumentation the 
fractional part of the carrier phase may be measured with an internal precision of 1 to 
3 mm, and with an overall accuracy (including the influence of the atmosphere) at the 
1- to 4-cm level, it is quite difficult to locate accurately the exact cycle of the carrier 
whose phase is being measured, to the effect that by an incorrect determination of the 
integer cycle ambiguity, an equivalent range error of a multiple k of the wavelength of 

the carrier may result; i.e., k x 19 cm for the Li frequency (or k x 25 cm for the L 2 
frequency). 

The estimation of this indeterminacy is usually carried out in the stage of the 
adjustment of the available observations by introducing an additional unknown bias to be 
estimated along with other nuisance parameters. The remaining problem in this 
approach stems from the difficulty in pinpointing the exact integer value (as the number 
of cycles can only be an integer) of the so-determined real-valued ambiguities. So far, 
many different approaches are used to resolve this ambiguity problem (see Langley et 
al., 1984; Bender and Larden, 1985; Bock et al., 1986), but as yet there is no efficient 
one for adjusting real and integer numbers simultaneously, not to mention the 
complications arising from the aliasing of the cycle ambiguities with the effects of 
multipath, atmospheric delays and orbital errors, particularly for data collected over 
long baselines. In this latter case, one can only allow these ambiguity parameters to 
assume real values. 

(b) The ambiguity problem is further compounded by the frequent occurrence of 
cycle slips. This is a commonplace phenomenon observed in GPS data when a satellite 
signal is obstructed in any way and can no longer be tracked. When signal lock is 
resumed, the fractional part of the measured carrier beat phase would still be the same 
as if tracking had been maintained all along. The integer number of cycles, however, 
exhibits a discontinuity, or cycle slip, which effectively adds a new ambiguity to the 
data. Such occurrences can be quite frequent, including many instances of loss of lock 
not recognized by the receiver (e.g., by the quality and tracking mode parameters of its 
measurements report), such as is shown, for instance, in Fig. 2.6 for the data of a 
typical observing session with TI-4100 receivers during the March '85 NASA/JPL 
experiment. 

Therefore, a vital question related to the preprocessing of GPS data is how well is 
one able to determine the limits of reconstruction of the number of whole cycles when 
phase lock is lost and regained; and hence, how well can one recover from such cycle 
slips by using various models for the actual change in range due to satellite motion and 
the phase error introduced by imperfect clocks. This initial assessment of the data is 
particularly important as it eliminates any damaged data which would otherwise alias 
the results of the orbit determination and baselines. 


2.3.1 Detection and Removal 


There several possible approaches to dealing with the cycle slip problem. An 
often-used approach (Hilla, 1986) is to perform a positioning adjustment in which the 
station positions are held fixed, and to edit the data by inspecting the corresponding 
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residual series for abrupt discontinuities. This has been shown to work, but it can only 
recover large cycle slips (less than a few tens of cycles). A similar approach is to use 
between-receivers, -satellites, and -time triple-difference observables (which are not 
affected by receiver or satellite clock errors) to determine station locations and 
subsequently, once convergence of station coordinates has been achieved, to search 
through the triple-difference residuals to identify large discontinuities in the double- 
difference observables which form the specific triple differences (Remondi, 1985). 
This can be tedious and time consuming. A common approach is to fit polynomials to 
various types of phase data (e.g. raw phase, phase with orbital motion removed, 
between-receivers phase, between-receivers and -satellites phase, etc.), and then use 
these polynomials to project the phase into the future. Using this approach Clynch and 
Coco (1986) have found, for instance, that at least a third-order polynomial is needed to 
model raw phase at the 1-3 cycle level for about 100 sec, whereas a second-order 
polynomial models at the same level about 180 sec of data from which orbit dynamics 
have been removed. These results varied as expected when orbit-removed data sets for 
various clock types were examined. Using observations spread over a 20-minute period, 
double-differenced carrier phases were predicted to within one half a cycle for up to 3 
minutes using crystal oscillators, 6.5 minutes using rubidium clocks, and 9 minutes 
using hydrogen maser clocks. 

In our analysis we used a similar approach based on a piecewise polynomial 
fitting procedure applied on data series of raw phase or between-stations single- 
difference measurements. Essentially, one may consider a time series of observed 
carrier phase measurements as composed of constituents that are of interest (e.g., 
satellite dynamic motion and satellite-to-receiver distance changes), and those noise or 
nuisance constituents (among them the cycle slips) which obscure the constituents one 
would like to study. The latter can be either periodic, rendering the series coloured, or 
some other type such as constant shifts or trends rendering the series non-stationary, 
or both. For example, from the discussion in this section it is easily understood that: 

- the occurrence of cycle slips observed in GPS data induces a datum shift in the 
observed time series which remains constant as long as carrier phase lock is maintained. 
That is, when signal lock is resumed, the fractional part of the measured carrier phase 
would still be the same as if tracking had been maintained all along. The integer number 
of cycles, however, exhibits a discontinuity or cycle slip which appears in the data as an 
abrupt jump in the carrier phase history (cf. Fig. 2.6). 

Touching for a moment on the aforementioned other aliasing effects on these cycle 

slips: 


- carrier phase measurements like all GPS measurements being intimately 
connected to precise timing, are also prone to the deterministic drifts of the satellite and 
receiver clocks involved in the measuring process. For the satellite rubidium and 
cesium atomic frequency standards, the latter are usually assumed to be linear or 
quadratic in time. Similar behaviour is often assumed valid for describing the long-term 
stability of the high-quality quartz crystal clocks of currently available GPS receivers: 

- a similar type of drift, which has also been identified in experimental 
observations, can be due to temperature variations inside the receiver (ibid); 

- multipath and imaging delays, which as already mentioned are usually noted 
when operating in a highly reflective environment, can also corrupt measured carrier 
phases at any elevation angle inducing periodic systematic noise in the observed time 
series. 
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All these non-stationarity-inducing constituents can be considered as 
systematic noise in the preprocessing stage, as long as we know their general form; that 
is, for example, the time at which cycle slips occur and the kind of trend (linear, 
quadratic, etc.)- The main point here is that an ideal preprocessing algorithm 
attempting to recover cycle slips in the GPS data should take into account the 
interrelationship of the hidden couplings exerted by these effects. Knowing as much 
about the constituents of the observed data series of carrier phase measurements, in the 
preprocessing stage, it is possible to use a piecewise continuous function of degree n 

n 

F(t) = 2 ajOj(t m )» <|> k (t m ) (2.18) 

j=1 

to approximate m carrier phases <|) k (or single differences) observed by a station (or 
two stations respectively if single differences are used) to satellite k where aj are 
coefficients to be determined and 


Oj s ( 0 O , 0! O n , cos G)jt, sincOjt ] (2.19) 

are base functions providing for datum shift biases, linear trends, quadratic (or other 
type) trends and periodic components, (Oj, of the kind described above. 

In practice, if this approach is to be applied solely for cycle detection and 
recovery purposes, only the coefficients pertaining to the datum shift base functions 
(0 O ) are of interest. If there are q cycle slips, present in the data series of a particular 
session, one ends up dividing the total session in sub-blocks Iq free of cycle slips, so that 

the size of cycle slips between blocks is determined by a series of a continuous function 
of the form (2.18) with identical coefficients except for the zero-order term; i.e., 


q n 

F(t) = 2 a oi + 2 aj 0j(t m ) 

i=i j=i 


( 2 . 20 ) 


so that an additional bias term is added to the original function for each sub-block I q . 

Although seemingly straightforward, the foregoing approach tacitly assumes that 
we know exactly the time at which the cycle slips occur in the data, that is to say, the 
division of the observation session into cycle slip free sub-blocks is known a priori. 
Such breaks, however, are often not obvious at the time of observation and can only be 
detected by a post-mission close examination of the data. In an effort to make this 
detection as automatic as possible, we have further implemented an algorithm which 
takes advantage of the short-term characteristics of the observed carrier phases to 
detect the occurrence of cycle slips. This is based on the fact that the change in GPS range 
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is exceptionally smooth (because of the extreme distance-3 to 4 earth radii-to the GPS 
satellites) and hence, highly predictable. Clock and propagation effects do not normally 
destroy this smoothness, except perhaps in the case of sudden ionospheric disturbances 
and highly dynamic receiver motion. 

Consider a sequence of carrier phase measurements <|>j, j=1,2 N. Let P 1 be 

the value at time t of the unique polynomial of degree zero (i.e. a constant) passing 

through the point (ti,<j>i); so Pi = 4>i (ti )- Likewise define P2, P3 Pn- Also let 

P12 be the value of the unique polynomial of first degree passing through both (ti ,<|) -j ) 

and (t2,<|>2)- Likewise define P23. P34. .... P(N-i)N- Similarly, one can define the 
higher order polynomials, up to P123..N' which is the value of the unique interpolating 
polynomial passing through all N points. This process can be schematically visualized as 
a tableau of various P’s with "ancestors" on the left leading to a single "dependent” at the 
extreme right; e.g., for N=4 


ti : 

<J>1 = p 1 

p i 2 



t 2 : 

<t>2 = P 2 


P 123 




CO 

C\J 

a. 

P 1 234 

(2.21) 

t 3 : 

<t>3 = p 3 

P 34 

P 234 


t 4 : 

a 

it 





Each unique polynomial P can be constructed by using the recursive relationship 
between a "dependent" and its two "ancestors" 


( M i+j) p i(i+1)...(i+j-1) + Or 1 ) p (i+1)(i+2)...(i+j) 

p i(i+1)...(i+j)= (2.22) 

(ti-t i+j ) 


which holds because the two ancestors already agree at points tj + i_ tj + j + i. A further 
improvement on this recurrence relationship is achieved by keeping track of the small 
differences 


c j,i - p i...(i+j) • p i...(i+j-1) 



(2.23a) 

D j,i = p i...(i+j) ' p (i+1)...(i+j) , 

j=1,2, .. 

., N-1. 

(2.23b) 


These differences satisfy the relation 
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Pj+1,i * Dj+1,i - C j>i+ i - Dj j 


(2.24) 


which is useful in proving that the recursive relationship in (2.22) can be converted to 
one involving only the small differences C and D, so that 


(tj+j+1 -t)(Cj ii+ i-Dj (i ) 

Dj + 1,i (2.25a) 

(V^i+j + 1 ) 

(tj-f) (Cj t i + i-Dj j) 

«V,J ■ < 2 - 25b > 

O.-t.+j+t ) 


In practice, at each level j, after each column in the tableau (2.21) is completed 
based on N given carrier phase values, if a predicted carrier phase value, 9 (t), at epoch 
t is sought, one can compute such a value by forking up or down the tableau which is 
equivalent to an interpolating or extrapolating value from the given <|>j values plus a set 
of C and D corrections which form the "straightmost" route through the tableau so that 
the partial approximations are centered insofar as possible on the target epoch t. The use 
of C and D corrections effectively makes the interpolation or extrapolation one order 
higher. The algorithm is easily implemented to follow the L 1 and L 2 signals and project 

the phase into the future. If the difference between the observed and projected values 
exceeds a pre-defined threshold, the algorithm interprets it as the occurrence of a cycle 
slip, so that an additional term is added in (2.20). The level of tolerance of these 
projections is obviously different for various types of clocks and the sampling rates of 
the data. From tests with TI-4100 data it has been found that 5 consecutive data points 
at 30-sec intervals were usually needed usually to model the raw phase at the 1 -3 cycle 
level for about 150 sec, which compares well with similar findings by Clynch and Coco 
(1986), who used algebraic polynomials and data from which orbit dynamics had been 
removed. This is considered adequate for the purpose of identifying the time of 
occurrence of suspected cycle slips in the data. 

The foregoing process was applied to several data sets prior to the analysis 
leading to the results to be reported in Section 5. In each case a system of equations 
(2.20) was solved simultaneously in a least-squares procedure similar to the least- 
squares spectral analysis described by Wells et al (1985) to estimate the coefficients 

a 0 j and aj. For the base functions Oj, Chebychev polynomials were found to perform 
better than algebraic polynomials simply because the Chebychev approximation is very 
nearly the same polynomial, which among other polynomials of the same degree, yields 

the smallest maximum deviation from the true function <}>(t). The degree of the 
polynomials required for cycle slip detection depends mainly on the length of the data 
series and the sampling of the data, as well as the smoothness (or noisiness) of the data 
to be approximated. Experience, for instance, with several data sets from the Spring 
’85 campaign has demonstrated that a degree 6 or 7 was typically required in order to 
produce piecewise fits of the raw phase over a data span of 1 hour with an rms error at 
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the 1 -cycle level or better, whereas a degree of 3 or 4 was adequate when fhe satellite 
dynamics were removed prior to the fit. These results varied as expected when data sets 
collected with different clock types were examined. Fig. 2.7 is a typical example of the 
fits obtained from the single-difference data from the Owens-Mojave baseline on day 90 
when both stations were operating with cesium clocks. In this case, no periodic terms 
were necessarily included in the fitting process to recover the cycle slips present in the 
original data series shown in Fig. 2.6. 


2.4 Clonk Error s 


In GPS, like VLBI, the tracking data being recorded are based on one-way signal 
travel time and hence, inherent to the measurements are the errors due to the receiver 
clock and those due to the clocks of the transmitting satellites. Although the rubidium 
and cesium frequency standards used onboard the GPS satellites are considered more than 
adequate for most applications of GPS, the requirements of geodynamic applications 
demand that biases in both the satellite and receiver clocks be properly modelled or 
otherwise be accounted for. As explained at the beginning of this section, a common way 
to eliminate or greatly reduce these biases is through double differencing or equivalently 
by including separate clock bias parameters for each satellite and receiver clock at each 
measurement epoch (which is equivalent to a "white noise" process). The differencing 
techniques are usually applied to a set of "quasi-simultaneous" measurements from 
different stations and satellites made over such a short time interval that the errors due 
to the clocks can be considered constant. In reality, however, one cannot disregard the 
fact that the estimation problem at hand (be it orbit determination or baseline 
estimation, or both) involves measurements over longer intervals of time, which in 
turn is forcing us to make some assumptions about the behaviour of the clock errors as a 
function of time. For example, in order to increase the degrees of freedom in the 
adjustment of the observations, a bias, a drift, and a drift rate are sometimes used (i.e., 
a quadratic clock model is assumed). This might be a satisfactory model to use for the 
satellite clocks and it is easy to handle, but is not generally accurate for describing the 
behaviour of the receiver clocks (except when a short observing period is involved or a 
high-precision frequency standard like a hydrogen maser is used). For the most precise 
applications, it is generally necessary to try to eliminate these errors altogether 
through direct modelling as opposed to allowing them to be absorbed by other 
parameters. 

One possible way of doing so, is to use a continuous time-state representation 
from which it is possible to derive (albeit at the cost of some added computational effort) 
the corresponding discrete time-state representation at the observation times. From 
experience gained so far with GPS operations, it is generally accepted that the frequency 
standards used in the currently available GPS geodetic receivers are quite stable; 
however, their behaviour is not deterministic, but rather undergoing stochastic 
variations in both the time and the frequency domains. From empirical verification of 
experimental work (e.g., Heroux, 1986), it has been determined by examining 
residuals after fitting a quadratic polynomial to GPS interstation clock errors that the 
residuals from such fits have random walk characteristics which are actual changes in 
the frequency behaviour of the clocks; hence, such variations can be interpreted as a 
random walk in frequency (which gets integrated into time), and an independent random 
walk in time. Stated differently, the variances of these two random walks can always be 
associated with the general clock characteristics such as long-term frequency stability 
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and short-term time stability. This makes a continuous model more suitable than any 
discrete time model since GPS measurements in most receivers are sampled at intervals 
from 6 sec to 5 min during each observing session, while the random walk behaviour is 
essentially continuous. The practical implication of this sampling process is that at the 
observation epochs there is a correlation between the random input to time and the 
random input to frequency caused by the continuous integration of the random walk in 
frequency. This is often disregarded for the sake of more simplified models. 

The form of the continuous process model underlined by such characteristics can 
be described as (e.g., Gelb, 1974; Jones and Tryon, 1987) 


dX(t) = F X(t) dt + D dt + dw(t) (2.26) 


with 


X(t) = 


X(t) 

p - 

0 1 

W(t) 

* r — 

0 0 


(2.27) 


and 


dW(t) 


e (t) 
T1 (t) 


o 

d 


(2.28) 


where x(t) is the clock error state, w(t) denotes the frequency random walk component, 
d denotes the deterministic frequency drift and e(t) and T)(t) are independent white noise 
processes with spectral densities c e 2 and o^ 2 respectively. For a discrete time interval 

8t, the prediction of the state vector X(t+8t) is 


where 


<D(8t) s I + F 8t 


= <2>(8t) X(t) + 

I 8 ' 

1 8t -x 
0 1 

J 0 

- 0>(5t) X(t) + 

D(5t) 

_ 1 1 

8t 1 



d T 


(2.29) 


0 1 


and 


D(5t)= d 


* 2 
8t / 2 

8t 


(2.30) 


(2.31) 
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with the covariance matrix of the error introduced by the random input to this 
prediction given as 


Q(5t) = J 


f 8t 

1 8t -x 


a 2 o 
e u 


1 0 

* 0 

0 1 


o ° 2 

v ^ 


8t -x 1 


dx (2.32a) 


2 3 2 
8t a £ + V 9 8t 


2 2 

8t 


2 2 

^8t 


8t a 




(2.32b) 


which illustrates how, although the white noise processes involved are independent, 
prediction over even a finite time interval introduces a correlation between time and 
frequency through Q. 

In a GPS data reduction scheme there are usually several clocks involved. The 
foregoing formalism allows the state equation for N clocks to be obtained by extending the 
two-variable state vector of each clock into a hypervector of 2N elements, and the 
corresponding transition matrix O to a 2N x 2N matrix consisting of 2x2 blocks (of the 
form (2.30)) on the main diagonal and zero off-diagonal blocks. Similarly, the 
covariance matrix Q of the random input error is extended to a 2N x 2N block-diagonal 
matrix with the 2x2 blocks on the diagonal being of the form (2.32a or 2.32b) where 
a e 2 and o^ 2 generally may be chosen to be different for each clock. Going one step 

further, during the adjustment of the observations, usually one clock is used as the 
reference clock and the differences between each clock in the observing session and the 
reference clock are inferred from the single-difference observations. This differencing 
produces an observation vector y(t) 


y(t) = A(t) X(t) + v (t) (2.33) 


where A(t) is the design matrix of the measurement partials with respect to the clock 
state vector and is of the form 
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1 

1 


0 

0 


(2.34) 


A(t) 


0 - 1 
0 0 

0 0 


0 0 

0 - 1 

0 0 


0 ... 

0 ... 

0 ... 


0 
0 

- 1 0 


and has dimensions (N-1) x 2N, and v(t) represents a random vector of observational 
errors assumed to be normally distributed with zero mean and an associated covariance 
matrix R independent of the input random noise dW and assumed to be diagonal with 
constant elements equal to the variances of the observational errors (usually known 
from the electronics of the receiver). 
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SECTION 3 


ESTIMATION OF GPS ORBITAL PARAMETERS 


3.1 General Considerations 


In general, every observation t to a GPS satellite (be it pseudorange, integrated 
Doppler or carrier beat phase) is a function of the satellite’s position at the observation 
time; i.e., 


t = C(z, p, x) + noise 


(3.1) 


where x is a vector of parameters of interest such as station coordinates, time delays, 
etc., p is a vector of parameters which describe the forces acting on the satellite, and z 
is the vector of orbital parameters which describe the motion of the satellite around the 
earth in a unique way. However, the choice and the number of such parameters are by no 
means unique. 

According to Newton’s second law, the satellite motion is defined as a particular 
solution of a system of second-order differential equations: 


r" « «F ( t ; r.r'.p) 


(3.2) 


where r = r(t) is the position vector of the satellite in an inertial Cartesian frame, r ' 
is the satellite’s velocity vector, and r" is the acceleration vector. The parameters p 
describe, for instance, the gravitational attractions of the earth, sun and the moon, solar 
radiation pressure, etc. Their choice and their number depend mainly on the lengths of 
the satellite arcs considered and as such are not unique. Equation (3.2) does not define 
the satellite’s orbit uniquely either. Fundamentally, the procedure for determining the 
orbit of the satellite can be considered as the process of assuming the satellite to be 
under the influence of a known force field and then using observations such as t to 

determine which solution to the equations of motion one should "choose". This is because, 
even if the forces acting on the satellite were known, an infinite number of solutions to 
the differential equations in (3.2) exist until boundary conditions are imposed-such as 
values for the six components of the initial position and velocity, six Keplerian 
elements, or similar elements at some chosen epoch. 

Let z„ be a vector of such a set of initial conditions for the satellite state and p 0 
and x c be vectors of initial estimates for the parameters p, and x respectively. 
Generally, the orbit estimation model consists of two sets of functions • | and £F 2 



& i (r.p.x ; i) = 0 
&2( r >*o) = 0 


(3.3a) 


(3.3b) 


where £F-| represents the purely geometric relationship between satellite positions and 

observables, whereas £F 2 represents the orbital motion model relating the satellite 

positions at any epoch and the initial state vector. If the initial estimates used are 
reasonably "good"; i.e., the components of the force model used to define the equations of 
motion are sufficient to yield a computed orbit close to the "true" orbit, their 
corrections would be proportionately small so that it may be possible to obtain them by 
solving a first-order approximation of the non-linear estimation model (3.1) based on 
linearized observation equations of the form 


St - @£/az) z=Zo 8z + 0£/3p) p . Po Sp + (dt/dx) x=XQ 8x (3.4a) 

= D z t Sz + D p t Sp + D x t Sx (3.4b) 

with D z , Dp, D x denoting the partial derivatives with respect to x, p and z respectively. 
In the sequel the following notation will be used extensively to denote: 

DyX ... the partial derivative of a scalar x (or vector x) with 

respect to a scalar y (or vector y) ; 

d y x ... the total derivative of a scalar x (or vector x) with 

respect to a scalar y (or vector y) . 

Evidently, any variations of the initial state vector errors 8z 0 would give rise to 
corresponding variations in the satellite state errors Sz related by 


8z = D Zq z 8 Zq + D p z 8p 


(3.5) 


where the components of D Zo z and D p z are solutions of a special form of the linearized 

equations of motion known as the variational equations and represent the sensitivities of 
the orbit to perturbations in the parameters z 0 and p 0 . Consequently, Eq. (3.4) takes 

the form 


St = D z t D Zq z 8z 0 + D z t D p z Sp +D p t Sp + D x t Sx . (3.4c) 
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As will be discussed in some detail in Section 5, there can be variations in the way these 
parameters are treated in an adjustment process which, in turn, define different orbit 
solution strategies or modes of analysis, each with different capabilities to reduce 
particular systematic errors and improve the orbit and baseline determinations. A 
completely general parameter estimation process should be able to solve for the best 
estimates of any combination of the orbital parameters 

z oi1 ’ z oi2’ z oi3’ z oi4, z oi5’ z oi6, >2 n s 

Pi > P2> •••> Pn> 

where n s is the number of satellite arcs being observed, and n is the number of dynamic 

parameters of the force model. Such complete generality is not always necessary 
however. If for instance, only relatively short arcs (up to one revolution or 12 hours 
duration for the case of the GPS satellites) were to be processed, it would allow 
modelling of W in equation (3.2) using relatively few parameters, p n , and would solve 
or update the set of state vector parameters z 0 j. Moreover, one would assume that some 

or all of these parameters (e.g., low-degree potential coefficients, lunar and solar 

gravitations, etc.) are partially known a priori, a fact that can be easily implemented 
through appropriate variance/covariance information in the adjustment process. Often, 
even this set may contain more "free" parameters than required. If the observations 
originate from a relatively small area and if the number of stations observing 
simultaneously is small, it will probably not make sense to solve for all z 0 \ elements 

either; instead, one would perhaps only be able to solve for one element (e.g., 

responsible for a possible along-track error) per satellite arc with any degree of 

certainty. 


3.2 Equations of Motion 


In the satellite equations of motion (3.2) r and r* can be given in any convenient 
system. In this section the equations of motion will be derived in a geocentric, inertial, 
Cartesian coordinate system as defined in detail in Section 4. Neglecting relativistic 
effects, the Newtonian equations of motion in the inertial frame are given as 


r"(t) = V V(r,t) + a(r,r',t) 


(3.6) 


where V h D r is the gradient operator, a is the vector of all disturbing accelerations 

acting on the satellite and W is the gradient of the earth’s gravitational potential 
referenced in the inertial frame; hence, its explicit functional relationship with time so 
that 


V ( r ,t) = V(r) + Vi(r,t) 


(3.7a) 
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where V(r) is the gravitational potential referring to an earth-fixed frame moving 
relative to the true inertial frame. Expressing V(r,t) in (3.6) by the way of Eq. (3.7a) 
results in a convenient decomposition of the gravitational potential in a radial part 

V(r) = (p/r) V V(r) = - (pr/r 3 ) (3.7b) 

where p is the product of the gravitational constant G and the mass of the earth M, and a 
small disturbing part Vi(r,t) for which 

VV^r.t) « V V(r) (3.8) 

so that a regrouping of the terms a(r,r',t) and W^r.t) as 

p(r,r',t) = a (r.r'.t) + VV^r.t) (3.9) 

leads to a homogeneous and a perturbed decomposition part of Eq. (3.6); i.e., 

r"(t) = V V(r) + p(r,r\t) (3.10) 

which essentially describes the central problem of dynamic satellite geodesy. 


3.2.1 The Homogeneous or Two-Body Problem 


The solution of (3.10) where p(r,r',t) s 0 is known as the homogeneous or the 
two-body problem and leads to the well-known result that the trajectory of the motion 
is a conic where the particular form (i.e., ellipse, parabola or hyperbola) depends on 
the initial conditions. Since the GPS satellites are in elliptical or near-circular orbits, 
the geometric interpretation of the solution is illustrated in Fig. 3.1. The osculating 
ellipse is the two-body orbit the satellite follows on the orbital plane with the center of 
mass of the earth as one of its focii. Since the gravitational force field is radially 
symmetric in this case (i.e. r and r’ are collinear), the motion is confined to the plane 
which is tangent to the true orbit at any instant of time t. The stationarity of the orbital 
plane, the perigee, and the size and shape of the orbit together with the constant period of 
the motion lead to a realization that the Keplerian motion is fully describable in terms of 
six parameters, of which one is a function of time. The choice of these parameters is not 
unique. A particular set of parameters, the familiar Keplerian elements, are 
particularly common for the study of the satellite perturbations. They are defined as 

shown in Fig. 3.1. The right ascension (Q) and inclination (i) define the orientation of 
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Fig. 3.1 - (a) above, the satellite orbit as projected onto a unit sphere, 
(b) below, the geometry of the orbital ellipse. 


35 







the orbital plane in space; the argument of perigee (co) defines the location of the 
perigee in the orbit; the semi-major axis (a) and eccentricity (e) define the size and 
shape of the orbit respectively, whereas one of the anomalies is normally used to define 
the position of the satellite in the orbital ellipse; i.e., the true anomaly (f), or the 
eccentric anomaly (E) which relates to the true anomaly by 


tan(f/2) = [(1+e)/(1-e)]^ tan(E/2) (3.11) 


and the mean anomaly M, which would be characterized as the angular distance from 
perigee if the satellite were moving at a constant velocity equal to the mean motion n 

n = (n a -3 )^ = d t M(t). (3.12) 

That is, 

M = n (t-Tp), Tp ... time of perigee passage (3.13a) 


usually given explicitly through Kepler's equation 


M = E - e sinE 


(3.13b) 


which can be solved via a series of expansion or iteration methods. 

These six elements, which may be considered as the six integration constants of 
the homogeneous equations of motion, are uniquely related to x,y,z and x'.y’.z' at any 
instant of time and can be used instead of them, the advantage of doing so being mainly 
that they do not vary nearly as much along the orbit (with the exception of the 
anomalies) as the Cartesian coordinates and velocities do. The transformation from the 
Keplerian representation to the position and velocity vectors is given as (Kaula, 1966) 


r = r(x,y,z) = r(a,e,i,co,Q,M) 

(3.14a) 

= R 3(-^) r i (■') Raf-w) q 

(3.14b) 

r'=r'(x,y,z) = r'(a,e,i,( 0 ,Q ,M) 

(3.15a) 

= R 3 (-Q) Ri(-i) R 3 (-co) q' 

(3.15b) 


where Rj(0) are the usual rotation matrices for a rotation angle 0 , and 
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q s q[qi.q2.<?3l T 


(3.16a) 


and 


with 


and 


q’ = q , [q’ 1 ,q , 2,q , 3] T 

(3.1 6b) 

qi = a (cosE - e) 

(3.17a) 

q 2 = a (1-e 2 )^ sinE 

(3.17b) 

O 

II 

CO 

cr 

(3.17c) 

q’l = - n a sinE / (1-e cosE) 

(3.18a) 

q '2 = n a (1-e 2 )^ cosE / (1-e cosE) 

(3.1 8b) 

q' 3 = o. 

(3.18c) 


3.2.2 The Inhomogeneous or Perturbed Problem 


In order to obtain a solution of the perturbed problem of the equations of motion, 
it is necessary to recast the non-linear vector second-order differential equations (3.2) 
to a non-linear system of three scalar differential equations of second order whereby the 
six integration constants are now considered as functions of time; i.e., 


r s r(x,y,z) s r(a(t),e(t),i(t),co (t),fl (t),M(t)) (3.19a) 

r’ sr’(x,y,z) = r'(a(t),e(t),i(t),co (t),fl (t),M(t)) (3.19b) 


This is the familiar method of variation of arbitrary constants which conveniently 
bypasses the problem posed in (3.2) that a direct solution in terms of Cartesian 
coordinates is intractable. Simply put, since the effect of small perturbations acting on 
the satellite will cause the satellite orbit to deviate slowly from an unperturbed 
Keplerian orbit, r and r' will also still correspond uniquely to an orbital ellipse that 
instantaneously corresponds to the perturbed motion, hence the problem reduces to 
solving for the time variations of the Keplerian elements instead. 
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Lagrange obtained this transformation of the equations of motion in an ingenious 
but rather complex way described in detail in Brouwer and Clemence (1961) and Kauia 
(1966). The end product is a system of six first-order differential equations 


6 

2 [ Sj» s k ] d ( s k * DgjR , j=1 , 2 ,.. .,6 (3.20) 

k=1 


where d t s k denotes the derivative of the s^ element with respect to time; R is the 

perturbing force function whose form will depend on the particular type of perturbing 
source (e.g., earth oblateness, lunar-solar attractions, etc.); s k or Sj represent the 

Keplerian elements a,e,i,co,Q,M, and [Sj,s k ] are the Lagrange brackets defined as 


6 

[Sj.Skl = 2 [ D Sj r D Sk r' - D Sj r' D SR r ]. (3.21) 

k=1 


A desirable property of the Lagrange brackets that facilitates their evaluation is their 
invariance with time , so that once evaluated at any convenient point in the orbit, such 
as, for instance, at the perigee where E=0, equation (3.20) can be inverted to yield the 
Lagrange planetary equations of the time derivatives of the Keplerian elements; i.e., 


d t a - (2/na) D M R (3.22a) 

d,e = [(1-e 2 )/(na 2 e)] D M R - [(1-e 2 )^ /(na 2 e)] D W R (3.22b) 

d t co = [- coti/(na 2 (1-e 2 ) ,/ * )] D,R + [(1-e)^/(na 2 e)] D e R (3.22c) 

d t i = (coti/(na 2 (1-e 2 )^)] D m R - [1/(na 2 (1-e)^ )sinij D n R (3.22d) 
d t Q - [1/(na 2 (1-e)'^sini)] D|R (3.22e) 

d t M = n-[(1-e 2 )/(na 2 e)] D e R - [2/(na)] D a R (3.22f) 

which are of the general form 

d t s k = s(a,e,i,co,Q,M; p, z 0 ) (3.23) 

where z c is the vector of the initial conditions at t 0 
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Oo)*^ (to).M(to)}. 


(3.24) 


z 


3.3 Likely Causes Shaping Orbital Errors 


3.3.1 Effects of Errors in the Initial Conditions 


As with all geodetic satellites the orbit determination problem for the GPS 
satellites in theory involves a complex process of integration of the equations of motion, 
but in practice ultimately reduces to an orbit improvement problem using observations 
which are themselves functions of many satellite positions, station positions and other 
"nuisance" parameters like satellite and receiver timing delays, atmospheric 
propagation delays and so on. As already mentioned, the formulation of the differential 
equations of motion as just outlined leads to an infinite number of solutions even if all 
disturbing forces are known; that is, until boundary conditions are imposed-such as the 
values for the initial position and velocity. The GPS observations are then used to 
determine as accurately as possible these initial conditions. Consequently, errors in the 
GPS orbits can arise from errors in the forces acting on the satellites and errors in the 
computed boundary conditions. In this context, the geopotential and a few other forces 
like the solar radiation are usually considered as the sole source of error in the satellite 
forces, and tracking station locations as the main source of error in obtaining erroneous 
boundary conditions. In principle, errors in the tracking station locations can be 
investigated separately from errors in the satellite forces, although in practice it is 
frequently difficult to completely separate orbit errors into those directly related to the 
satellite motion and those related to the station positions. This is because both are 
largely interrelated: to determine the orbits one needs to know both the geopotential and 
the positions of the tracking stations; and to obtain the position of other stations with 
observations, it is necessary to know the orbit and, thus, the forcing field. 

To understand the effect of an incorrect initial state on a computed orbit, it is 
necessary to look back at the perturbed equations of motion in (3.10) written as 

P 

r" + r-A (3.25) 

r 3 

where A denotes the disturbing acceleration caused by the perturbing forces acting on 
the satellite. When a first-order perturbation is applied on (3.25), we obtain 


P 3p 

8r" + 5r 5rr = 8A (3.26) 

r 3 r 4 

where 8A is considered to represent the perturbing acceleration resulting from either 
an incorrect set of geopotential coefficients, the initial state, or the force model. The net 
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effect on the computed orbit becomes quite clear if the resulting satellite position error 
vector 8r is expressed in the orthogonal coordinate system {S,T,W} of the radial, along- 
track and out-of-plane components shown in Fig. 3.1(a) that moves with the satellite in 
the inertial frame and has its origin on the satellite as defined by the reference orbit. If 
{ 82;, 8t|, 8£ } are the corresponding components of the error vector 8r in this frame, 
then 

8r= 8^? + 8q f) + 8C C (3.27) 

where £ , fj and £ are the unit vectors in the radial, along-track and out-of-plane 
directions, which notably are time dependent with 


d t $ * d t M i 

(3.28a) 

d t fj « -d t M T) 

(3.28b) 

dtC ■* 0 . 

(3.28c) 

In this frame 


8A = D r 8R t + (1/r) Du8R T) + (r sinu)' 1 Dj8R £ 

(3.29) 

where 8R represents errors in the force function, u-f+co is the argument of latitude, 
and i is the satellite inclination. Substitution of (3.27) into (3.26) yields in view of 
(3.28) 

2\i 

d 2 t 8£ - [ (d t M) 2 + ] 82; - 2 d t M d t Sri - d 2 t M 8ri 

r 3 


= D r 8R ■ 8S 

(3.30a) 

d 2 t 8r| +[ (d t M) 2 ] 8ti + 2 d,M d t 8£ - d 2 t M 8£ 

r 3 


= (l/r) D u 8R 5 8T 

(3.30b) 
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d 2 t 8C + 8£ - (r sinu)' 1 Dj8Rs8W (3.30c) 

r 3 

where 8S, 8T, and 8W are the transformed errors in the acceleration components defined 
in the {S,T,W} frame. Neglecting the small eccentricities of the GPS satellites, equations 
(3.30) can be further simplified as 


d 2 t 8£ - 3 n 0 2 8£ - 2 n 0 d t 8r| = 8S 

(3.31 a) 

d 2 t 8ri + 2 n 0 8^ - 8T 

(3.31b) 

d 2 ,8C + n 0 2 8£ = 8W 

(3.31c) 


known as Hill’s equations (cf. also Colombo, 1985) which interestingly reveal that 
across-track errors are independent of the radial and along-track ones. 

The general solution of these equations can always be considered as being 
composed of a particular solution of the inhomogeneous equations, that is, when the 
equations (3.31) include terms explicitly dependent on the gravitational potential or on 
other external forces such as solar radiation, and a general solution of the homogeneous 

part of the equations (i.e. 8R b 0) which reflects how orbital errors Aa(t 0 ), Ae(t 0 ), 

Ai(t 0 ), Aco(t 0 ), AQ(t 0 ), AM(t 0 ) in the initial state at t 0 will propagate as functions of 
time. The complete solution is obviously their sums. 

The form of the induced errors in the orbit from an incorrect initial state can be 
obtained by solving the linearized homogeneous equations of motion with the errors of 
the initial state as initial conditions, i.e. 


8a(t) - 

1 Aa(t 0 ) 

(3.32a) 

8e(t) - 

1 Ae(t 0 ) 

(3.32b) 

8i(t) * 

Ai(t 0 ) 

(3.32c) 

8co(t) . 

= Ao)(t 0 ) 

(3.32d) 

800) 

= AO(t 0 ) 

(3.32e) 

8M(t) 

= -(3/2) (n 0 /a) Aa(t 0 ) (t-t Q ) + AM(t 0 ) 

(3.32f) 


or equivalently, in the radial, along-track and across-track system 
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8£(t) = d t 8^(t 0 ) (sin n 0 t/n 0 ) 


- [(2/n 0 ) d t 8ri(t 0 ) + 3 8^(t 0 )] cos n 0 t 

+ [(2/n 0 ) d t 8Ti(t 0 ) + 4 8£(t 0 )] (3.33a) 

8ri(t) = (2/n 0 ) d t 8^(t 0 ) cos n 0 t 

+ [(4/n 0 ) d t 8r|(t 0 ) + 6 8£(t 0 )] sin n 0 t 
+ [5r| (t 0 ) - (2/n 0 ) d t 8^(t 0 )] 

- [ 3 d t Sri(t 0 ) + 6 n 0 8£(t 0 )] t (3.33b) 

8C(t) = 8C(t 0 ) cos n 0 t + d t 8C(t 0 ) (sin n 0 t/n 0 ) (3.33c) 


and consequently 

d t 8£(t) = d t 8^(t 0 ) cos n 0 t + [2 d t 8ri(t 0 ) 

+ 3 n 0 8£(t 0 )] sin n 0 t (3.34a) 

d t 8ri(t) = - 2 d t 8£(t 0 ) sin n 0 t 

+ [(4 d t 8ri(t 0 ) + 6 n 0 8£(t 0 )] cos n 0 t 

- [ 3 d t Sri(t 0 ) + 6 n 0 8S(t 0 )] (3.34b) 

d t 8C(t) = - n 0 8£(t 0 ) sin n 0 t + d t 8C(t 0 ) cos n 0 t . (3.34c) 


As to be expected, if there is an error in the period of the satellite motion, the 
satellite along-track error grows linearly with time, whereas the satellite altitude/ 
radius exhibits an error which will not average to zero and hence, will likely translate 
directly into the determination of the station positions in a systematic way. These terms, 
although they arise mathematically from a solution of the homogeneous perturbed 
equations of motion, should not be construed as trivial additions to the perturbed 
satellite motion from the physical point of view. By contrast they illustrate the very 
resonant character of the linearized equations of motion. Like the orbits of the altimetry 
satellites (also usually having eccentricities approaching zero) for which Colombo 
(1984a) showed that the same expressions also hold, the GPS orbits can behave like an 
undamped dynamic system whose unforced response produces oscillations in the satellite 
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motion at the system’s natural frequencies: 0 and n 0 . When driven by a disturbing 
force that has such a periodic component at the natural frequencies, the forced response 
of the system will include similar oscillatory terms which grow larger as their 
frequencies approach 0 or n 0 (cf. Fig. 3.2). 


3.3.2 Effect of Unmodelled or Improperly Modelled Forces 


The particular solution of (3.31) reflects the contribution of the forces acting on 
the satellite. Although 8A may well be thought of as being due to errors in the 
geopotential resulting from errors 5C nm and 8S nm in the harmonic coefficients, these 
are not expected to be significant for (the very well known to date) low degree and order 
terms which are usually required in the case of the GPS satellites for modelling the 

earth’s gravitational force during the orbit integration. Rather it is likely that 5 A 
would be due mainly to unmodelled or improperly modelled forces, particularly of the 
non-gravitational nature like the solar radiation pressure. In his more recent work 
Colombo (1985) related his earlier developments on this aspect from the altimetry and 
the satellite-to-satellite tracking satellites to the GPS case suggesting that the same type 
of forced response may occur under the influence of non-gravitational forces, such as: 

- the solar radiation effects whereby the solar radiation pressure resulting from 
the electromagnetic interaction between radiation from the sun and the surface of the 
spacecraft, may repeat itself at each revolution if the orientation of the spacecraft with 
respect to the sun is maintained (which happens anyway with the GPS satellites). Any 
departure from perfect periodicity will be due to the imperfections of the force model, 
the slow changes in the relative geometry of the earth, the sun, and the satellite, which 
would reasonably produce some resonant orbit error close to the fundamental frequency 
n 0 ; 


- the variable part of the earth’s gravitational field arising from tidal and other 
deformations of the solid earth and the oceans. This likewise depends on the relative 
geometric arrangement of the spacecraft, the sun or the moon and the earth, and so may 
cause resonance effects for similar reasons; and 

- the constant force along the spacecraft’s y axis caused by the lack of perfect 
alignment of the solar panels (Fliegel et al, 1985) known as the y-bias effect. This is 
causing periodic changes (given the way the satellite orientation is maintained in space) 
almost at a rate of once per revolution which is also likely to excite resonance. 

If errors in the initial state and any unmodelled (or incorrectly modelled) forces 
are large enough, their resonant effects will have a simple mathematical approximation 
expressing the complete solution of (3.31) in the form (Colombo, 1984b; 1985) 


K K J 

Auj(t) ■ £ A ki t k cosn 0 t +£ B k j t k sin n 0 t +X CjjtJ (3.35) 

k=0 k=0 j=0 
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where i=1,2,3 and Auj(t) s {A^.Arj.A^}. In fact, Colombo (ibid) has verified with 

simulations that with this parameterization, some terms will effectively be very close 
to zero, and others will be related in a simple way, thus leaving at the most, 18 
parameters or about 6 per component defining 


A£(t) = cos n 0 t + sin n 0 t + t cos n 0 t 

+ t sin n 0 t + E^ (3.36a) 


Ar)(t) = Aq cos n 0 t + Bq sin n Q t + Cq t cos n 0 t 

+ Dq t Sin n 0 t + Eq + Fq t + Gq t 2 (3.36b) 


AC(t) = cos n 0 t + B^ sin n 0 t + t cos n 0 t 

+ t sin n 0 t + E^ + t (3.36c) 


from which similar expressions for the velocity components can be easily derived by a 
further simple differentiation. Numerical results with actual data using this 
formulation in the present work will be reported in Section 5. 

Having outlined the adverse effects that may arise on the GPS orbits by the likely 
errors in the initial state and the force models, it is in order that the mechanisms which 
shape them now be considered. 


3.4 Perturbations to GPS Orbits 


The GPS satellites are designed to provide precise three-dimensional position, 
velocity and time information anywhere in the vicinity of the earth. Perturbations due 
to: 


- the gravitational attraction of the earth; 

- the gravitational perturbations of the moon, the sun and the planets (i.e. the 
third-body effects); 

- solar radiation pressure; 

- atmospheric drag; 

- magnetic forces; and 

- the variable part of the earth’s gravitational potential arising from tidal and 
other deformations of the solid earth and the oceans, 

act to drive the satellites away from their operational orbits. Because of the high 
demand for accuracy, the Ground Control takes care that the satellite subtracks are 
maintained over the lifetimes of the satellites. This orbit maintenance is taking place in 
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the form of small satellite manoeuvres to periodically correct the orbit and maintain the 
repetitive nature of the satellite ground tracks. 

The major perturbations on the GPS satellites act on the spacecraft in a variety 
of ways and timescales, not all of which are important or relevant to the determination 
of the orbits usually used for the determination of geodetic baselines; these are usually 
of up to 1 week in duration. Typically, the disturbing forces just outlined may be 
considered according to the characteristic amount of time over which they act; i.e. 
generally, 

- secular perturbations which are constant over time (e.g. gravitational 
perturbations due to the earth’s oblateness); 

- short-period perturbations superimposed on a precessing Keplerian 
ellipse which act on periods of a few hours; i.e., the orbital period of the 
satellites (e.g. gravitational perturbations for the low-degree harmonics of 
the geopotential); 

- long-period perturbations which act over periods of a few days to 
several weeks (e.g. the third-body effects); and 

- resonant perturbations arising from the repetitive ground-track 
condition as a result of particular combinations of the tesseral and sectorial 
harmonics with the given secular rates in o>, Q, and M; and the forced 
response effects induced under the repetitive influences of such causes as those 
outlined in Section 3.3.2. 

As it turns out, the orbital elements of the GPS satellites undergoing variations of 
interest are co, Q, M and a. Of these, the three angular ones undergo secular, long- 
period and resonant perturbations. The semi-major axis, a, does not vary secularly or 
in long periodic terms, but is influenced almost solely by resonant perturbations which 
change its value by a few meters per day at the most, but in turn, is the primary cause 
of a nodal longitude acceleration analogous to the one experienced by geosynchronous 
satellites. 


3.4.1 Disturbing Function of the Gravitational Potential 


3. 4.1.1 Linear perturbations 


The integration of the Lagrange planetary equations as given in (3.22) requires 
that the partial derivatives of the disturbing function R be evaluated with respect to the 
Keplerian elements. The development of the disturbing function due to the geopotential 
begins with the standard representation of the earth’s gravitational potential in 
spherical harmonics C nm , S nm (up to a maximum degree N) 


V « V„ m (r, ♦•,».) 
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N n 

= (GM/r) {1 + I I (a E /r) n P nm ( cos<|>’) 
n=2 m=0 

[C nm cos mX. + S nm sin mX.]} (3.37) 

in terms of the the mean earth radius a E , the geocentric radial distance r, the geocentric 

latitude <f>' and longitude X. of the satellite and the associated Legendre functions P nm 

(Heiskanen and Moritz, 1967). For the treatment that follows, equation (3.37) has to 
be expressed in terms of the Keplerian elements and 0, the Greenwich sidereal angle that 
accounts for the earth’s rotation, in the form (Kaula, 1966) 

R(a,e,i,o),Q,M,t) = I I V nm (a,e,i,a),fl,M,t,0) (3.38) 

n m 

where 

n °o 

V nm (a f e,i,co,£2,M,t,6) = (GM/a n+1 ) I F nmp (i) I G npq (e) x 

p=2 q=0 

^ Spnipq(®t^i»M,0(t)) (3.39) 

with 

/ C nm> n-m even 

Snmpq = cos[(n-2p)co + (n-2p+q)M + m(fl-0(t))] 

\-S nm , n-m odd 

/ S nmi n-m even 

+ sin[(n-2p)co + (n-2p + q)M + m(£l-0(t))] (3.40) 

\ C nm , n-m odd 

The functions F nmp (i) and G npq (e) are known as the inclination and eccentricity 

functions respectively, which involve the evaluation of binomial coefficients and 
factorials. Ordinarily these functions (and their derivatives) can be very troublesome 
and computationally expensive to evaluate especially for high degrees, n, of the harmonic 
expansion of the potential. For the near-circular GPS orbits two important properties of 
G np q(e) somewhat alleviate this problem: that they are approximately proportional to 

etal* so that they converge rapidly for small eccentricities and with increasing q; and that 
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1 if q « 0 


/ 

Lim G npq (e)« 

e -» 0 \ 0 otherwise 


(3.41) 


This is a tremendous simplification since the inner summation of (3.39) is restricted to 
a single term; namely, q-0. Furthermore for the GPS satellites, because they are in 
high stable orbits at about 20,000 km altitude they are much less affected by the short 
wavelength features of the geopotential. Thus, only a subset of the spherical harmonic 
coefficients (usually up to degree and order 8) is used in the integration process which 
simplifies the computation of the eccentricity and inclination functions with any of the 
many analytical expressions found, for example, in Allan (1965), Kaula (1966), the 
recurrence relations in Kostelecky (1985) or the Fast Fourier Transform relations 
proposed by Goad (1987). Here the notation in the disturbing function in (3.38) by the 
way of (3.39) and (3.40) has been chosen to correspond closely with that of Kaula 
(1966). Thus G n pq(e) is identical to the expressions given in Kaula, whereas the 

F nm p(i) can be evaluated from the expression given by Allan (1965) as 


(n+m)l 

Pnmp(0 «= 

2 n p! (n-p)l 


2n-2p 2p 

X (-l)k ( )( ) c 3n-m-2p-2k s m-n+2p+2k (3.42) 

k k n-m-k 


(c = cos 14 i and s s sin 14 i) which is somewhat more simple. 

The representation in (3.39) of the potential admits a rather elegant form for 
solving the Lagrange planetary equations when the angular rates of a, e and i are 

constant, while co, Q and M are linear functions of time. Careful inspection of (3.38) 
indicates that the single largest effect in R by the gravitational field would be produced 
by the earth’s oblateness; i.e., the second zonal harmonic term C 20 , causing most of the 
secular variations in the mean Keplerian orbit. Consequently, by introducing 


R - (p/a) C 20 (ae/a) F 201 (i) G 210 (i) 


(3.43) 


for the disturbing function into the Lagrange planetary equations, the secular variations 
of the Keplerian elements would be written as 


d t a - 0 


(3.44a) 
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d t e = 0 (3.44b) 

d t i = 0 (3.44c) 

d t co - y A w C 20 [(a£/a) 2 / (1-e 2 ) 2 ] (1-5 cos 2 i) (3.44d) 

d t Q = 3n C 20 [(a E /a) 2 / 2(1-e 2 ) 2 ] cos i (3.44e) 

d t M = n {1- 3C 20 [(a E /a) 2 / 4(1-e 2 ) 3/2 ] (3cos 2 i-1)}. (3.44f) 

To first order the change in s k s {a.e.i.co.Q.M} over a complete orbit is 

Ais k = ^ d t s^ dt , (3.45) 

with the integration taking place along the reference orbit where a, e and i are assumed 
constant during the time of integration and the radial distance of the satellite is the same 
as for an unperturbed orbit, and 

(0 = co 0 + d t co (t-t 0 ) (3.46a) 

Q = Q 0 + d t Q (t-t 0 ) (3.46b) 

M = M 0 + d t M (t-t 0 ) (3.46c) 

9 = (Q O -0 O ) + (d t Q-C0 E ) (t-t 0 ) (3.46d) 

where t 0 is the time of the initial conditions and oo E is the rotation rate of the earth. The 
final result is: 

2^/(p/a)(a E /a) n F nm p(i) Gnpq(®)( n ' 2 P +c l)^nmpq(®) 

^1 a nmpq = (3.47a) 

[(n-2p)d(C0 + (n-2p+q)dtM + m(d t Q-co E )] 
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V(n/a 3 ) (a E /a) n f*nmp(i) ®npq(®) ^ 

X (1- e 2) ! / M(1-e 2 ) !/ Mn-2p + q)-(n-2p)]S nmpq (©) 

A 1 e nnpq= 

e[(n-2p)dtC0 + (n-2p+q)dtM + m(dtfi-co E )] 

V(n/a 3)(a E /a) n F nmp (i) x 

X G npq( e )K n " 2 P) cosi * m ) S nmpq( 0 ) 

^l'nmpq = 

(1-e 2 )^ sini[(n' 2 p)dtC 0 + (n-2p+q)d|M + m(dti2-co E )] 

V(^/a3) (a E /a) n [(1-e 2 )^ e-iF nmp (i) d e G npq (e) - 
-coti (1-e 2 )'^ djF nmp (i) G npq (e)] S nmpq (0) 

^1 ^nrnpq = 

[(n-2p)d t o) + (n-2p+q)dtM + m(dt£2-ci) E )] 

V(p/a 3 )(a E /a) n djF nmp (i) G npq (e) J S nmpq (@) 

e 

^nmpq 

(1-e 2 )^ sini [(n-2p)d t co + (n-2p+q)d t M + m(d t Q-co E )] 


(3.47b) 


(3.47c) 


(3.47d) 


(3.47e) 
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Fig. 3.2 - Perturbation effects on the GPS satellite orbits 
due to the second zonal harmonic of the geopotential. 


50 


AlMnmpq 


(3.47f) 


V(p./a®) (ag/a) n Fnmp(’) [-(1-e^) e ^ d 0 G n pq(e) - 

-(2n+1) G n pq(e)] J* 

Snmpq (0) 

e 


[(n-2p)d t co + (n-2p+q)d t M + m(d t Q-C 0 E)] 


where 0 * [(n-2p)co + (n-2p+q)M + m(Q-coE)], and djF nmp (i) and d e G npq (e) are the 
derivatives of the eccentricity and inclination functions with regard to i and e 
respectively. 

The typical perturbations (3.47) are short-period effects impressed upon a 
processing Keplerian orbit. The combined secular and periodic perturbations are quite 
significant even for short orbital arcs. These vary from about 800 m in i to about 5000 

m in Q. When transformed into the {S,T,W} coordinate system attached to the satellite, 
which is closer to giving a direct measure of the satellite position error, these effects 
are as shown in Fig. 3.2 for a 2-day arc. The perturbations due to the higher harmonics 
of the gravitational potential are found to be quite significant, but are attenuated quite 
rapidly because of the high altitude of the GPS orbits. Fig. 3.3 shows, again resolved into 
radial, along- and across-track components, the growth of the effect of the higher-order 
terms up to degree and order 4 over the 2-day arc. The resulting combined satellite 
errors exhibit a similar time dependence to the errors described by Eq. (3.34). This 
illustrates that over orbital arcs it is indeed possible to "soak up" most of the satellite 
position error due to this geopotential error by appropriate adjustment of the satellite's 
initial orbit parameters. 


3.4.1. 2 Higher order (non-linear) perturbations 


From the equations already derived for the mean rates of change over one 
revolution, one can obtain the effect of the first- and second-order gravitational 
perturbations by integrating A-jSk with respect to the given element over its cycle 
period. Approximating the second-order solution by a Taylor’s series expansion 


dtSk = dtSk(So) + D S jSk AiSj + ... 

■ d t (A! s k ) + d t (A 2 s k ) (3.48) 
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Fig. 3.3 - Perturbation effects on GPS orbits due to higher 
degree and order terms of the geopotential excluding C 20 


in which, upon subtracting out the first approximation d^A^k), yields 


d t (A 2 s k ) = D Sj s k AiSj (3.49) 

from which the second-order perturbations are computed as 

A 2 s k = I X D Sj s k A-|Sj ] dt . (3.50) 

t 


Normally, finding a second-order solution by evaluating (3.50) would be 
computationally inconvenient since the evaluation of A 2 s k would require a considerable 

amount of differentiations and multiplications. Fortunately, again because of the near- 
circular GPS orbits, matters are greatly simplified. This is because, as it turns out 
from a careful inspection of the terms involved in (3.50), only the changes in the mean 
anomaly (arising from the perturbations in the semi-major axis), the node and the 
perigee are significant-in addition to the movement of the earth; furthermore, the 
first-order perturbations which should be included in (3.50) require the effects of only 
the terms R 2 oio °f the disturbing potential for each A-|Sj which entails considerable 
simplification. 
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Fig. 3.4 - Schematic representation of the earth's 
ellipticity causing resonant effects on the GPS orbits. 


3.4.2 Resonant Perturbations 


The possible existence of resonance is well known in geodetic orbit theory. 
Resonance occurs for orbits whose periods have some simple ratio to the length of a day. 
For such orbits the satellite is influenced by the same disturbing gravitational feature at 
repeated intervals and resonance is thereby induced in the satellite motion. The GPS 
satellites experience a type of resonance analogous to the type experienced by 
geosynchronous satellites arising from the tesserai and sectorial harmonics, 
particularly C 22 and S 22 - This is illustrated in Fig. 3.4. In a reference system rotating 
with the earth, it is clear from symmetry that when the satellite is off either principal 
axis of the earth's equatorial ellipse, the bulge nearer the satellite has a stronger 
attraction on the satellite due to its shorter distance from it. If, for instance, the bulge 
is "ahead" of the satellite, it tends to exert a pull along its orbit thus transferring 
angular momentum to the satellite which, in turn, should result in an increase in the 
period and potential energy, and a decrease in the mean linear and angular velocities. 
Because of the conservation of angular momentum, the increase in the latter results in 
an increase in the semi-major axis of the orbit. A similar type of resonance (but one 
harder to visualize in such simple terms) exists for the GPS satellites arising from 
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particular combinations of tesseral and sectorial harmonics. Because of the high altitude 
of the orbits the GPS satellites tend to have low angular velocities, the satellites slowly 
drift in time by some 4 minutes a day which appears as an eastward shift of the satellite 

subtracks by approximately 1°. Because the commensurability of the GPS satellite 
orbital periods with the earth’s rotation rate is not exact, the longitude of the ascending 
node has a secular rate upon which resonant effects are impressed. 

Generally, the corresponding variational equation for this nodal regression can 
be derived from the equation of the earth-referred longitude of the node crossing given 
generally for any satellite by 


X,g = (1/Sq)(co + M) + ft - 0 


(3.51) 


where s 0 is the ratio of two integers in the relationship s = s 0 + As where s is the 
node-to-node satellite revolutions per earth rotation and As<1. For the GPS satellites 
obviously s 0 = 2, and As * 0.012. Note X s is given as the sum of angles lying in two 

different planes and should not be confused with the conventional geodetic longitude 
measured along the equator. From (3.51 ) 


X s = ^ [(co + M) + ft - cog] 


(3.52) 


or 


d(X s ■ K [(djG) + djM)] + djft - co E 

- (n/2) [1-(2/n 2 a)]D a R - !4 [(1 -e 2 )/(n 2 a e)] D e R + 


+ {(1-^cosi)/[n 2 a(1-e 2 )^]} DjR - co E (3.53) 

in view of (3.22). 

The form of the disturbing function for the resonant effects can be derived from 
equation (3.38) by considering the limits on the summations to be restricted only to 
terms in which 


!4 [(d t co + d t M)] + d t ft - co E = 0; (3.54) 

i.e., restricting the indices n, m, p and q to those satisfing the conditions: 
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(i) 


n - 2p = m/2 ^ [n-(m/2)] even , 

p = (2n - m) / 4 


(3.55a) 


(ii) q = 0 (3.55b) 

due to the small eccentricity of the GPS orbits, so that such terms as V 3210 . v 44io> 

• • • • 

V 5220 etc - will give rise to terms with an integral multiple of \'/t ( 0 ) + M) + Q - co E ]. 

Because of the presence of the (a E /r) n term in the disturbing function (3.37), the 

dominant effect will rise from the harmonics for which n is small. A further 
simplification is possible by making the observation that 

^nm = (1/m) arctan (S nm /C nm ) (3.56) 

so that the disturbing function for the resonant perturbations becomes 

R “ ^nm = (H/ a ) £ (a E /a) n V(C nm + S nm ) X 
n,m 

X Fnm,n-(m/2)(e) cos[m(co + M + Q - 0 - X nm )]. (3.57) 

Because of the presence of the mean anomaly M in the disturbing function, there will be 
a major perturbation of the semi-major axis from the first of the equations (3.22), i.e. 

d t a = - (2/na) I (a e /a)" V(C nm + S nm ) x 

n,m 

X Fnm,n-(m/2)(e) m sin[m(0) + M + Cl - 0 - X nm )]. (3.58) 

During integration of the Lagrange planetary equations, 

d t \y = m d t X s - q d t to (3.59) 

(denoting y = m (X s - X nm ) - q (0) is explicitly assumed constant by considering that 
the only variations are the secular rates d t M, d t co and d t Q due to oblateness and the 
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rotation rate of the earth. The resultant perturbations of the Keplerian elements are 
sinusoidal with periods that are submultiples of the orbital period, earth rotation, and 
any combination of these. Consideration of Eq. (3.55) however, reveals that d t y must 
also be constant for any term V nmpq that satisfies Eq. (3.55). This implies that other 

resonant terms must be added to those already defined. These correspond to those with 
periods P c (usually termed "circulation" periods) related to the orbital period P 

through the resonant parameter R c 


R C = P C /P=«1 +(s 0 /As) (3.60) 

- 1 + (2 / As) 

recalling that s 0 = 2 for the GPS satellites. Generally, resonance can be observed 

physically for all satellites making more than 10 revolutions per day, whereas for those 
satellites for which P < 1 0 revolutions per day, it is possible to observe resonance only 
when the orbit period is nearly commensurate with the earth's rate of rotation (i.e. As 
is small); put differently, when the resonant parameter R c is less than 10, resonant 

effects can be ignored. For the GPS satellites given that As = 0.012, R c reaches values 
in the hundreds (R c = 161), so that the resonant perturbations begin to appear as forced 
along-track oscillations given as 


A = a (AM + Aco + A Q cosi) 


(3.61) 


impressed on the motion of the secularly processing ellipse. AM, Aco and AQ in (3.61) 
are evaluated with Eq. (3.47) simplified by the assumptions (3.41) directly applicable 
to near-circular orbits, and fixing m to m=s 0 for the fundamental resonance or m = k s 0 

for other resonances. 

In addition, resonant perturbations would occur when the divisor in (3.47) 
approaches zero, i.e. 


(n-2p)d t co + (n-2p+q)d t M + m (d t £i - ©e) = 0 . (3.62) 


A particular case relevant to the GPS satellites falling into this class occurs near the 
critical inclination i=63.4° for which 


© « 0 . 
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Currently the prototype Block I satellites have inclinations in the vicinity of the critical 
inclination so perigee motion is very slow. Block II satellites are planned to be near 

55° inclination, so that the perigee advance should be at a faster rate. 


3.4.3 Lunisolar (Third Body) Perturbations 


The gravitational perturbations caused by a third body on the satellite's orbit are 
treated by defining a function R which is the third-body potential. This takes the form 


R d = (GM m d /r d ) { [1 - (2r/r d ) cosy + (r2/r d )]-'^ (r/r d ) cos V } (3.63) 


where m d and r d is the mass and the geocentric distance of the third body, y is the angle 
between the geocentric vectors r to the satellite, and r d to the third body respectively. 
For the case of the GPS satellites, only the sun and moon effects are usually treated in 
this context. 

Generally for short satellite GPS arcs (<i/3 revolution or = 4 hours duration), 
one can safely assume that the third body is stationary during that period. In this way, 
treating the earth, the satellite and the third body as point masses, the gravitational 
perturbations resulting from the moon and the sun can be developed analytically in a 
manner very similar to the case of the gravitational perturbations. The form of the 
disturbing function in this case takes the form (Kaula, 1962 ; Buffett, 1985) 


n max n n 

R = (G m d /r d ) I (a/r d ) n III Fnmp( e ) x 
n=2 m=0 p=0 q=-» 

X Hnqp(®) S nmpq^’^'M) (3.64) 

where 


/ A nm> n-m even 

S*nmpq = cos[(n-2p)(0 d + (n-2p+q) M d + m Q d ] 

\-B nm , n -modd 


/ B nm> n-m even 

+ sin[(n-2p)co d + (n-2p+q) M d + m Q d ] . (3.65) 

\ A nm> n-m odd 
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Here H npq (e) are the second inclination functions (in place of G npq (e) for the earth 
gravitational case), and A nm and B nm are spherical harmonics given as 


A nm - [(n-m)l/(n+m)l] e m P nm ( sin8 d) cos ma d 

(3.66a) 

B n m * [(n-m)l/(n+m)l] e m P n m(sinS d ) sin ma d 

(3.66b) 

/ 1 , m=0 


e m ■ 

(3.66c) 

\ 2, m?t0 



with a d , 8 d and {a d ,e d ,i d ,co d ,£2 d ,M d } being the declination and right ascension and the 
orbit elements of the third body. The general expressions for the H nqp (e) functions can 
be found, for example, in Buffett (1985); they need not be repeated here. 

For a more complete treatment of the third-body disturbing function, this study 
has followed Buffett’s (1985) suggestion incorporating the motion of of the third body 
following an earlier development by Giacaglia (1973) whereby the spherical harmonics 
A nm and B nm are rotated in the ecliptic reference frame resulting in Keplerian elements 

of the third body which are referred to in the ecliptic frame and hence can be more 
precisely approximated using the mean elements of the IAU 1980 Nutation Theory 
(Wahr, 1981; Seidelmann, 1982). In this way, the difficulty in approximating the 
equatorial elements of the sun and moon as simple functions in time is eliminated. The 
final result is a very complex-looking disturbing function 

n max n n n n oo oo 

R - (G m d /r d ) I (a/r d ) n I I I I I I 

^nmp(') 

n= 2 m»0 s--n p-0 t*0 q»-<*» u«-« 

^ Pnst('d) ^npq(®) 

* U nms Snmsptqu^’^’M) (3.67) 


where 


U nms - ^nms , S £0 (3.68a) 

U*nrns = M) n ' s [(n-s)1/(n+s)!] U nms , s < 0 (3.68b) 
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r 2 n+m n-m 

U nms = I (-1) n-m-r ( m+s+r )( r ) cos 2r+m+s (e/2) X 


r ” r 1 


x sin 2n ‘ 2r ‘ m * s (e/2) 


(3.68c) 


with 


r 1 = max {0;-m-s} 
r 2 = min {n-s;n-m}. 


and 


S nmsptqu ■ 


r !4 { cos [ {(n-2p)to+(n-2p+q)M+mQ} + 

{(n-2t)© d +(n-2t+u)M d +s(fl ( j+JC/2)}] + 

cos [ {(n-2p)o)+(n-2p+q)M+mQ} - 

{(n-2t)d) d +(n-2t+u)M d +s(Q d +JC/2)}] } ... s even 

(3.69) 

Vz { sin [ {(n-2p)co+(n-2p+q)M+mQ} + 

{(n-2t)0) d +(n-2t+u)M d +s(Q d +it/2)}l + 


sin [ {(n-2p)co+(n-2p+q)M+mQ} - 
{(n-2t)© d +(n-2t+u)M d +s(Q d +Ji/2)}] } 


... s odd 


for the zonal harmonics, and 


f cos [ {(n-2p)fl)+(n-2p+q)M+m(Q + Jt/2)} - 


{(n-2t)0) d +(n-2t+u)M d +s(Q d+ JC/2)}] ... (m-s) even 


Snmsptqu 


(3.70) 


(-1) n * s sin[ {(n-2p)to+(n-2p+q)M+m(Q + Jt/2)} - 


L {(n-2t)« d +(n-2t+u)M d +s(Q d+ ii/2)}] ... (m-s) odd 

for the non-zonal harmonics. 
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The magnitude of the lunar perturbations even for short GPS orbital arcs can be 
of the order of 50-500 m. By comparison, given that the value of m s / r s 3 for the sun 

is » 0.46, the perturbations of the solar third-body effects are about half of those of the 
moon. The gravitational attractions of the planets are even much smaller than those of 
the sun, so that their effects on GPS satellites are usually considered negligible. The 
indirect effect on GPS satellite motion caused by tidal deformations of the solid earth and 
the oceans due to the sun and the moon, is insignificant for short orbital arcs. For longer 
arcs (> 2 days duration), perturbing accelerations due to the ocean tide of the order of 

10‘ 9 m/s 2 can excite resonant effects resulting in orbital perturbations of over 2 
meters. These can be accounted for by way of empirical accelerations described in 
Section 3.3.2. 


3.6 Solar Radiation Pressure 


Radiation from the sun carries with it not only energy but also momentum. When 
this radiation is absorbed by some object like a satellite, the momentum transferred 
from the radiation creates a repulsive pressure upon the impact of the light protons on 
the satellite. Generally, the pressure normal to the incident direction is 


P s «(I/c) 


(3.71) 


where I, the intensity (solar flux) of the solar radiation at the mean earth-distance 
from the sun (i.e., one astronomical unit), is approximately 1360 Joule/m 2 sec and c is 
the speed of light in free space, so that P s = 4.6X1 O' 6 N/m 2 . The resulting force on the 
spacecraft's surface is proportional to the effective area A of the surface normal to the 
incident radiation 


F r= C R A 


(3.72) 


where Cr is the surface reflectivity constant. The acceleration acting on a spacecraft of 
mass m at a distance r ss ( r ss = r - R sun ) from the sun is 


— 2 
h r r sun 

r B sp rss (3.73a) 

m r3 ss 


CrA 

= P S 

m 



'ss 


3 


(3.73b) 


60 



where r sun (=149 598 500 km) denotes an astronomic unit and F r is directed opposite 
to the satellite-sun line. 

The GPS satellites have irregular, complex surface geometries and are 
constructed of materials with some surface portions reflecting diffusely and other 
portions reflecting specularly. In addition, they experience thermal changes in the 
surface area of the satellite that is illuminated by the sun as the satellite orbits the 
earth. Hence, in reality, at any moment the satellite will experience a net force that is 
not necessarily along the satellite-sun line plus a torque. A practical implication 
arising from such a complex satellite shape is that it is rather inappropriate to assign a 
single Cp value (especially for long arcs) and A for the GPS satellites. Clearly, a solar 

radiation model which was to take into account the separate contributions of the 
satellite’s main body, the solar panels, the antenna array and the engine assembly would 
require a very complex model like the ROCK-IV model (Flieger et at, 1 985). For short 
arcs a simpler "flat-plate model" with a single Cr estimated from observations would be 

quite adequate. 

A practical solution for the effects of the solar radiation pressure requires 
several simplifying assumptions. The discontinuity of the force due to the shadowing 
effects of the earth as the satellite enters and exits the shadow zone (Fig. 3.5) is taken 

into account by the inclusion of a shadow function factor v so that 


Cr A r ss 

r sp “ v Ps r ^sun . (3.73c) 

m r3 ss 


The development of an analytical theory in the presence of a shadow function is quite 
complicated for the reason that points of discontinuity are functions of the orbital 
elements and of the sun’s position. Kaula (1966) has treated the problem in an 
iterative manner using a simple case of a cylindrical shadow model. Lala (1971) 
developed an elegant semi-analytical model designed for differential orbit improvement 
programs utilizing high-accuracy observations, and as such, is particularly suitable for 
the case of the GPS satellites as our implementations have shown. The shadow function 
used is of the form 


v =!4[1 + sinx (1-cos 2 x)'^ ] , X “ (3.74) 


where X is the satellite angular distance from the center of the earth’s shadow and <E> is 
the angular distance between the shadow center and the shadow boundary intersection 
with the satellite orbit, with all angles being measured at the earth’s center; 
consequently, when the satellite is in the sunlight portion of the orbit 


X s [O.Jt] v =1 


(3.75a) 


61 



Fig. 3.5 - Solar radiation geometry 

and when the satellite is in the shadow portion of the orbit 

X€[-jc/ 2,0]-» v -0. (3.75b) 


Unlike the gravitational field case where the usual method of obtaining the 
change in the Keplerian elements of the orbit with time was to first find the disturbing 
function and then use Lagrange's planetary equations, this approach cannot be used in the 
present problem since the force cannot, in general, be written as the gradient of the 
potential. Hence no disturbing function exists. Instead the force must be dealt with 
directly. Once the force is found, it can be used in the Gaussian form of Lagrange’s 
equations in the manner indicated below. 

The solar radiation pressure accelerations can be expressed in terms of 
orthogonal components (cf. Fig. 3.1a) ss [S,T,W] T to be used in Gauss’ form of 
Lagrange’s planetary equations. These components are derived from the Cartesian 
components r" sp of the satellite acceleration vector via a rotation: 


s = [S,T,W] t = R sr r” sp 


(3.76) 


where 
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Table 3.1 - Solar radiation pressure parameters for use in Gaussian form of 
Lagrange’s planetary equations (after Lala, 1971). 
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F = (a 2 /p.) r" sp 

(1) M e = max power of the eccentricity 


retained in the series 

Ai = P| 1 element of R rs 


A 2 = r 2 i element of R rs 

(2) XCj = - e XC M i>2 

Bi = W cosCO 

XSj = - e XSj_i ^ i>1 

B 2 = W since 
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Rsr - R 3 «0(t)) Ri(i(t)) R 3 (Q(t)) Ri(-e) R 3 (-X. sun (t)) 


= {r U } 


i,j = 1,2,3 


(3.77) 


with A. sun being the ecliptic longitude of the sun. 

For short satellite arcs, it makes sense to assume that the orientation of the orbit 
with respect to the sun is constant. In this way, the elements of the rotation matrix R sr 

can be treated as constant which in turn allows the components of s to be treated as 
functions only of the true anomaly of the satellite. Introducing (3.77) into (3.76) 
followed by the insertion of the latter into the Gaussian form of the Lagrange planetary 
equations leads to 


I M 

D t z = 2n FA { X XCj cos'f + X XSj cos'f sinf } , z s [a,e,i,co,Q,M] (3.78) 
i=0 i=0 


where n is the mean satellite motion and the parameters FA, XCj and XSj are given in 

Table 3.1 for each Keplerian element. The trigonometric terms cos'f are usually 
replaced by equivalent expressions in terms of the mean anomaly M using the Hansen 

coefficients X n m * k ; that is 


cos'f = X 0 °.' + X [X u °.' + X. u 0 *' ] cosuM . (3.79) 

u=1 

Under the foregoing assumption that the orientation of the orbit is constant with respect 
to the sun within the time interval of integration, the solution of the equations of motion 
would comprise only short-period and secular terms. Over longer arcs, the coefficients 
FA, XCj and XSj could not be considered to be constant with the integration, resulting 
instead in long-period perturbations. This effect would be best observed in the semi- 
major axis which does not undergo a long-period motion from any other source. This 
perturbation is (cf. Murphy and Felsentreger, 1966) 


El 

Aa sr = - 2(a 2 /p) P s (A/m) S cosE + T V(1-e 2 ) sinE 

E 2 


(3.80) 


where Ei and E 2 are the eccentric anomalies of the satellite at the exit from and the 
entrance into the earth’s shadow. For longer arcs, the assumption of constant orbit 
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orientation with respect to the sun can still be tolerated if the formulation in (3.36) is 
adopted to account for any unmodelled effects. 

For GPS satellites the acceleration of a satellite due to the solar radiation 

pressure is of the order of 10* 7 m/sec 2 ; for short satellite arcs the perturbation effects 
on the orbit can be of the order of 2-10 m, whereas they may reach values from less 

than 100 m to over 800 m after 2 days. 

The indirect solar radiation arriving at the satellite after being reflected from 
the earth’s surface; i.e. the albedo effect, may be considerable especially on the orbits of 
satellites with high surface-to-mass ratios like the GPS spacecraft. The effect of the 
albedo pressure on the GPS satellite orbits is considered to be of the order of 1-2% of 
the direct solar radiation effect (King et al., 1985), which incidentally is of the same 
order of magnitude as the effect of the earth and ocean tides, and hence, can be neglected 
for most applications or considered in later developments as required. 
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SECTION 4 


TIME, COOBDINRTE FRAMES, RND OBSERUABLES 


For the processing of GPS data we usually require results in three dimensions 
although the GPS system is inherently four-dimensional. In doing so, there are a 
number of coordinate and time frames which are inherent in these computations and 
hence, deserve particular attention. To begin with, there are various coordinate systems 
which are naturally involved in the determination of the satellite orbits. The satellites 
are moving with respect to an inertial reference frame and hence, a fundamental 
requirement in the processing of GPS data is the definition of an inertial coordinate 
system for the formulation and the solution of the equations of motion. Other coordinate 
systems are required in order to develop the perturbing forces acting on the satellites. 
The station coordinates on the other hand are usually expressed in an earth-centered, 
earth-fixed (ECEF) reference frame which is moving with respect to the inertial 
reference frame because of the earth’s motion. The measurements made by the GPS 
receivers are essentially timing differences of the time of transmission (in the satellite 
time scale) and the time of arrival (in the receiver time scale) of a particular 
electromagnetic signal transmitted by the satellite and therefore are affected by the 
timing errors in the clocks involved. Such measurements necessarily involve detailed 
and careful consideration of a number of time scales. In particular, time enters into GPS 
data processing in three different ways: through the time tags of the measurements, 
through the time tags of the satellite ephemerides, and through the measurements 
themselves. These are intimately connected with the definition of the observables and 
for that reason, this section deals with the brief description of the celestial and 
terrestrial reference systems used in the orbit computations. The time systems inherent 
in the GPS system and the GPS observables will also be discussed. 


4.1 Time and GPS 


In order to fully appreciate the role of time in GPS data processing, it is 
necessary to define the various time systems involved and their attendant time scales. 
Some of these definitions are standard. Others are particular to the GPS system and the 
GPS data processing. These will assist in the subsequent development of the observation 
equations and the formulation of strategies relating to the elimination of GPS clocks and 
bias elimination errors. 

In general, there are three time systems that are used in GPS computations: 
dynamical time, atomic time, and sidereal time. Like any orbital body the motions of the 
GPS satellites can be expressed most readily as a function of a dynamical time system as 
defined by the theories of Newtonian or relativistic mechanics. This is the fundamental 
uniform time system that defines the independent variable in the satellite's equations of 
motion. When we generate a set of four-dimensional coordinates of the GPS satellites- 
"ephemerides"-we implicitly use dynamical time. In this system, the motions of the 
satellites and the perturbing forces acting on them can be used themselves to measure 


66 



dynamical time. For all practical purposes, the dynamical time is equivalent to an 
Atomic Time System (Moritz and Mueller, 1986). The latter is more easily realized by 
the atomic clocks, e.g. the international network of clocks coordinated by the Bureau 
International de I'Heure-BIH, that maintain the International Atomic Time scale (TAI). 
The atomic time scale defined by such clocks is not itself a true or perfect time scale but 
is the most accurate measure of "true" time and is used to define the origin and the rate 
of other time scales such as (1) the Terrestrial Dynamical Time , TDT (which 
represents the dynamical time for motion of bodies within the earth’s gravitational 
field) being related to TAI by 


TDT = TAI + 32.184 sec ; 


(4.1) 


(2) the Universal Coordinated Time, UTC which runs at the same rate as TAI but is 
incremented by 1 -second jumps (leap seconds) when necessary to account for the 
slowing down of the earth’s rotation with respect to the sun and is being transformed to 
TAI by adding a specific integral number of seconds (23 at present), and (3) the GPS 
System Time. The GPS System Time used by the Control Segment is derived from UTC as 
maintained by the atomic clocks of the U.S. Naval Observatory. GPS time was set to UTC 
on 6 January, 1 980 but because the leap seconds that are introduced in the UTC time are 
not accounted for in the GPS time, there presently exists an integer-second difference 
between the two time systems. Currently this difference is 4 sec, i.e. 


GPS Time = UTC + 4 seconds 


(4.2) 


which implies a constant offset of 19 sec between GPS time and TAI, i.e. 


GPS Time = TAI - 19 seconds. 


(4.3) 


In addition to the atomic time scales which provide a uniform measure of time. 
Sidereal Time is required to relate the inertial and terrestrial systems. Universal Time, 
which is the most common form of sidereal time, is a rotational time scale affected by 
the irregularities of the earth’s rotation. By correcting observed universal time for 
the effects of polar motion, according to the definition of UT1, this time reflects the 
actual angular orientation of the earth with respect to the vernal equinox at that instant. 
UT1 is derived from observations made by BIH and is published as corrections to UTC, 
i.e. 


UT1 = UTC + DUT1. 


(4.4) 


Since UTC is related to GPS System Time, the transformation between UT1 and GPS time 
is accomplished by a simple shift in origin. 

The relationship between UT1 and the Greenwich sidereal time is defined as 
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GMST at O h UT1 = 24110 s .54841 + 86401 84S.81 2866 T u 

+ 0 s . 0931 04 T u 2 - 6 s . 2 x 10’6 T u 3 (4.5) 

where GMST is the mean sidereal time defined by the Greenwich Hour Angle of a 
fictitious sun Sy and the mean equator and equinox of date, and T u is the number of Julian 
centuries elapsed since JD 2451545.0 (36525 days of Universal Time elapsed since 
J2000=Jan. 1.5, 2000). 

Individual atomic clocks maintain their own time scale often referred to as clock 
time. This is true also for the clocks of the GPS satellites and receivers, each of which 
maintain a scale ( satellite clock time vs. receiver clock time) which is a representation 
of the "true" time (be it GPS time or UTC), but in practice it has some error with 
respect to it. The GPS satellites carry both rubidium and cesium frequency standards. 
The space environment is kind to the GPS clocks, but the latter are physically left to 
drift off the standard GPS time system. However, their performance is continuously 
monitored by the GPS Control Segment and the amount by which they have drifted is 
accurately known. For example, synchronization between the GPS satellite clocks is kept 
to within about 20 nsec by the broadcast clock corrections, and the GPS time scale is 
synchronized to UTC to within 100 nsec. Future experiments like the Advanced 
Clock/Range Experiment (ACRE) are designed to test compact H-Maser frequency 
standard capability for future GPS satellites. This, in addition to laser ranging capability 
at the cm-level (for this experiment Block I GPS satellite #12 is to be retrofitted with 
laser retroreflectors for 1992 launch (Buisson, 1987)) will allow simultaneous 
comparison of H-Maser, cesium and rubidium clocks in space and the capability to 
separate orbit position and velocity errors from satellite clock errors. 

As with any measurement system based on the motion of satellites such as SLR 
and LLR, NNSS/TRANSIT, as well as GPS, another time scale implicit in the 
measurements is one that relates to the time tags of the satellite ephemerides. This is 
often referred to as the Satellite Ephemeris Time or SET (not to be confused with the 
formal Ephemeris Time (ET) system defined by the orbital motion of the earth about the 
sun). Just as individual clocks can have differences with respect to a time scale such as 
GPS time or UTC, the motion of the satellite can be considered as defining a unique time 
scale which in general will represent the satellite ephemeris time but will have some 
errors with respect to it. Given, for instance, a set of satellite positions (say, from 
range observations) and an ephemeris that relates the satellite positions to SET, one can 
define a time scale close to the SET but different from it due to errors in the ephemeris 
and the satellite positions themselves. This realization of SET may be referred to as 
Satellite Orbital Time (SOT). SOT can be considered as an approximate measure of 
dynamical time and should not be confused with satellite clock time. 

The time system inherent in the GPS computations involves all the time scale 
relationships just outlined. The physical GPS observable, be it a pseudorange or a 
carrier beat phase, is essentially the difference between the reference phase of the 
satellite oscillator at the time of signal emission and the reference phase of the receiver 
oscillator at the time of signal reception. In rigorous terms this clock difference is 
composed of the accumulated sum of five components: the difference between station clock 
time and proper time (ideally UTC) at the station, the difference between proper time 
and coordinate time at the station (i.e. UTC-TDT), the difference between the coordinate 
time at signal emission and signal reception, the difference between coordinate and 
proper time at the satellite (TDT - UTC + At re |; At re | = relativistic clock correction), 
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and finally the difference between proper time and clock time at the satellite. In 
commonly used GPS terminology, these differences can be readily interpreted 
(disregarding for the moment the influence of atmospheric effects and instrumental 
delays): the first and fifth differences are nothing else but the familiar station and 
satellite clock errors respectively: the second and fourth differences are general 
relativistic time transformations; the third difference is the solution to the light-time 
equation connecting the signal emission and reception event; i.e., the geometric range 
between satellite and receiver defined by the epochs of these events. If during the 
modelling of the observations these relationships are rigorously taken into account, then 
the GPS computation time may be equivalent to any of the above time scales. For most 
applications certain time offsets between these time scales may be assumed negligible 
and ignored altogether, hence defining by default the computational time by the time 
scale of one reference clock (at the satellite or receiver level), by the average clock 
times of several receiver clocks or even (as is often done routinely in GPS data 
reductions) by the time offsets obtained from the GPS pseudorange navigation solutions. 

These concepts will hopefully became more apparent by seeing how these time 
scales are used in the development of the GPS carrier phase observation equations. 


4.2 Coordinate Frames 


Although the observables obtained from GPS receivers are instrumented 
differently from one receiver type to another, it can be shown that they are all, be it 
pseudorange, carrier phase or Doppler, functions of the instantaneous geometric ranges 
between the satellites and ground receivers and their time derivatives. This geometric 
range is by far the larger portion of the observed range, and also the main cause of 
complexity in building the observational model because of the numerous transformations 
which are necessary to relate the inertial reference frame used for the formulation and 
solution of the equations of motion and for relating the satellite position to the earth- 
fixed, earth-cantered reference frame which the station positions are referenced to. 

The coordinate system adopted for the solution of the equations of motion in the 
present work is based on the fundamental astronomic system defined by the geocentric, 
equatorial frame with the equator and equinox of J2000.0 and the 1976 IAU conventions 
and the 1980 Nutation Series (Leiske et al., 1977; Seidelmann, 1982). 

For the numerical solution, in order to evaluate the accelerations due to the 
earth’s gravitational potential, the inertial coordinates of the satellite are rotated into 
the average terrestrial system (i.e. the system in which the spherical harmonic 
coefficients are referred to) using a series of rotations 


r AT = S N P r i . 


(4.6) 


where S, N and P denote transformations to be explained below for polar motion and UT 1 , 
nutation and precession respectively. When the gravitational accelerations are evaluated 
in the AT system, they are transformed back into the inertial system using the transpose 
of the complete transformation matrix S N P. 
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In (4.6) the orientation of the AT system with respect to the true system of date 
is determined via S which represents the transformation 

S - R 2 (-x) R-|(-y) R 3 (GAST) (4.7) 


where x,y are the coordinates of the polar motion (readily obtained from published BIH 
values) and the Greenwich apparent sidereal time is related to GMST (Eq. (4.5)) by 

GAST « GMST + 8y cos(e + 8e) (4.8) 

where e is the obliquity of the mean equator of date given by 


C = 23° 26' 21.-448 - 46."8150 T 

- 5.-9 X 10' 4 T 2 + 1.-813 X T 3 (4.9) 

and 8e is the nutation in obliquity. 8e along with 8y, the nutation in longitude, and e (= 

e + 8e), the obliquity of the ecliptic, relate the mean and true systems of date, Fig. 4.1a. 
Following Wahr’s (1981) nutation theory for an elastic, elliptical and oceanless earth, 
8 e and 8y can be expressed as a series expansion of linear combinations of five 
fundamental arguments a;, i.e., 

106 5 

8e « Z Aj(T) cos [ Z kjj a,- } 
i=»1 j=1 

106 5 

8y « Z Bj(T) sin [ I kjj aj ] 
i-1 j-1 

where Aj(T), Bj(T), kjj are known coefficients, and the aj terms represent: 

- the mean anomaly of the moon 

a-) « l * 485866".733 + (1325 r + 715922".633) T 

-i- 31 ".310 T 2 + 0".064 T 3 ; (4.11a) 


(4.10a) 


(4.10b) 
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Fig. 4.1 - (a) above, Precession angles, 
(b) below, Nutation angles. 
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- the mean anomaly of the sun 


a 2 = t' = 1287099". 804 + (99 r + 1292581". 224) T 

- 0".577 T 2 + 0".012 T 3 ; (4.11b) 

- the mean argument of latitude of the moon (FsL-Q ) 


a 3 = F = 335778". 877 + (1342 r + 295263". 137) T 

+ 13". 257 T 2 + 0".011 T 3 ; (4.11c) 

-the mean elongation of the moon from the sun (Dsl_-L\ and L and L' are the 
astronomic longitude of the moon and the sun respectively) 

a 4 = D = 1072261 ".307 + (1236 r + 1 105601".328) T 

- 6".891 T 2 + 0".019 T 3 ; (4.1 Id) 

- the mean longitude of the ascending lunar node measured in the ecliptic 

a 5 = Q = 4501 60". 280 + (5 r + 482890".539) T 

- 7".455 T 2 + 0".008 T 3 ; (4.1 1 e) 

where 1 r = 360° s 1296000". An additional set of terms, although not as yet 
implemented in the present software, can be added to the nutation variations (4.10) to 
account for (a) the out-of-plane nutations (not originaly included in the IAU 1980 

Nutation series) which are identical to (4.10) except for the replacement of sin <-» cos 

106 5 

8e = X Cj(T) sin [ X kjj a; ] (4.12a) 

i=i j=i 

106 5 

8y = X Dj(T) cos [ X kjj aj ] (4.1 2b) 

i=1 j=1 
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and (b) for the free-core nutations with nominal period of 460 days (Sovers and 
Border, 1987). 

Applying these rotations defines the transformation between the mean system of 
date and the true system of date through the nutation matrix N 


N = FW-e) R 3 (-6y) FMe) . 


(4.13) 


The transformation from the mean equatorial system of date to the equatorial system of 
the reference epoch (e.g. J2000.0) is represented by the three equatorial arguments 

£ A , Z a and 0 A (Fig. 4.1b) 


C A = 2306" .2181 T + 0" .30188 T 2 + 0."017998 T 3 (4.14a) 

z A = 2306". 21 81 T + 1".09468 T 2 + 0."018203 T 3 (4.14b) 

0 A = 2004". 31 09 T - 0".42665 T 2 - 0."041833 T 3 (4.14c) 

defining the precession matrix P as 

P = R 3 (-90°-z a ) R !(0 a ) R 1 (90°-C A ) (4.15) 


which completes the model for the orientation of the earth in our adopted J2000.0 
system of reference. 


4.3 Obseruables 


GPS observations can always be thought of as some type of biased range. Code 
pseudorange or carrier beat phase observation equations contain various bias terms 
modelling the effects of the imperfections of the time scales used to generate and time tag 
these observables, atmospheric delays, etc. In the case of the carrier phase the 
corresponding observation equation contains an additional linear bias term-the initial 
phase ambiguity. A comprehensive summary of the simplified form of these equations 
and the various linear combinations has been given by Wells et al. (1986). A more 
detailed account of the mathematical development of these equations can be found, for 
instance, in Remondi (1984) and Wei (1986) just to mention a few. In this overview, 
only the carrier beat phase observation equation will be given for completeness, since 
this was the observable dealt with mostly in this study. 


73 



It is well known that the phase O of an ideal oscillator relates uniquely to the 

true time x (which provides the theoretical frame for deriving the observation 
equations) by 


x 

0>(t) = O(x 0 ) + j f dx 


(4.16a) 


where f is the time-dependent frequency of the oscillator. For an oscillator of high 
stability, over a short interval (x-x 0 ), (4.16a) can be written as 


O(x) = O(x 0 ) + f 0 (T-T 0 ) 


(4.16b) 


where f 0 is the nominal frequency of the oscillator. In practice, the carrier beat phase 

measurement at some time epoch T(x), measured at the receiver time scale, is based on 
the phase alignment of the receiver-generated signal with the incoming carrier signal 
without the knowledge of which cycle would represent perfect cycle synchronization (see 

Fig. 4.2). Hence, the total phase <\> consists of an integer number of cycles N(<|>;T 0 ,T) 
counted since the initial lock-on time T 0 , a fractional part Fr(<J)) and an unknown 
integer number N of cycles at the initial epoch T 0 , so that 


<t> - N(4>;T 0 ,D + Fr(<(» + N(T 0 ) 


(4.17a) 


= ^measured + N(T 0 ) . 


(4.17b) 


Ideally, if the receiver maintains continuous lock on the signal, N is a unique number for 
each satellite-receiver pair. Otherwise, every time the receiver loses lock, a cycle slip 
occurs which effectively introduces a new ambiguity term (see also Section 2). 
Generally the frequency of the receiver oscillator is not constant with respect to "true 
time" and thus the clock time will have a non-linear relationship with respect to it. For 
the nominal l-i (*1575.42 MHz) frequency a 0.6 nsec error in clock time will create 

approximately 1 cycle of phase error. Since phase is physically invariant, the received 
phase of the receiver-generated carrier signal at reception time T is the same as the 
phase of the satellite-generated carrier signal at transmission time t, and therefore 


7 4 


<t> = 4>(t) • <t>(T) 


(4.18) 
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Fig. 4.2 • Schematic representation of carrier 
beat phase observable generation. 


or after considering Eq. (4.16b) and the relationship of receipt and transmit times 


t + dt + (p/C) - T + dT + 8 atm , (4.19) 


the equation for the instantaneous carrier beat phase observable for one satellite, one 
receiver and one epoch is readily obtained as 


<|> = - (f/c) p(x) - f [dt(x)-dT (x)] - 8 atm + N (4.20) 

where any satellite and station frequency deviations from the nominal value have been 
absorbed in the terms dt and dT respectively; c is the speed of light in the vacuum; and 
8 a tm denotes the delay caused by the travel of the signal through the earth's atmosphere. 


75 



This is caused by a change in phase velocity but, as the effect is small compared to the 
geometric range, it can be modelled as a time delay, phase delay or simply as an increase 
in range (cf. Section 2); p denotes the geometric range to the satellite determined by the 
position vectors of the satellite and the receiver in the chosen coordinate frame. In 
(4.20) the time-dependent terms have been purposely expressed in terms of the "true 
time". This does not usually present any problems for the clock terms, if these are 

eliminated by differencing. The geometric range p, however, is largely dependent on 
using a correct time scale relating the time of transmission to the satellite ephemeris 
time from which the satellite position vector and hence (together with a priori 
information on the station’s position) the range is calculated. This calculation is usually 
done iteratively, with the accuracy of the calculation being dependent on the difference of 
the time of reception used (i.e. measurement time tag) and the time of reception 
expressed in terms of the ephemeris time. The effect of time tag error on the 
computation of the geometric range is readily seen in the Taylor series approximation 


p(T) = p(T+5t) = p(T) + 8 t p(T) + [ (5T)2 p(T)/2] (4.21 ) 


where p and p are the rate of change (Doppler) and range acceleration (negative of the 
Doppler slope). For the height of the GPS satellites, the rate of change reaches a 
maximum of about 800 m/sec, whereas the maximum value of the second derivative of 

the range is about 175 mm/sec 2 , or less than 1 Hz/sec. Clearly, if the tag error is less 
than one psec both of these terms can be ignored, otherwise additional terms will be 
required in the adjustment to account for this type of error: ideally one extra parameter 
for every receiver and every epoch. However, this obviously would lead to a 
prohibitively large number of parameters and hence is seldom done; instead, the 
technique of differencing between-receivers, -satellites and -time, and combinations 
thereof is often preferred (albeit at the expense of some loss of geometry, and 
redundancy of the data as well as the introduction of correlations). 

The combination of Eq. (4.16) through (4.21) leads to the complete form of the 
observation equation of a carrier beat phase measurement from station 1 to satellite j 


<(>ii - - (f/c) [pi(T i ) + St! pi(T! )j 


-f [dtiOM-dT^T^] - 8i 1iatm + NJ, . (4.22) 


The single-difference phase observable, defined as the difference of simultaneous 
measurements made by two stations, 1 and 2, to the same satellite j (with simultaneity 
customarily defined at the receiver time scale level), is readily obtained from (4.22) as 

<t> j 12 = <> j 2 ‘ <l> j 1 (4.23a) 
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= {- (f/c) lpj(T 2 ) + 8 t 2 pJ(T 2 )] 

- f [dti(tj 2 )-dT 2 (T 2 )] - 5i 2iatm + Nj 2 } 

(f/c) [pid^ + Sx, pirn)] 

- f ldti(ti 1 )-dT 1 (T 1 )] - 5i 1( atm + N*1 } (4.23b) 

• • 

= - (f/c) {[pi(T 2 ) - pj(T i )] + [5x 2 p j( T 2 ) - Sx, pj(T 1 )]} 
-f[dti(ti 2 ) - dti(ti-| )] + f [dT 2 (T 2 )-dT I (T t )] 

- [8i 2i atm -Sii.atml + 1^2’ N M- (4.23c) 


The satellite clock errors still appear in (4.23c) because the two observed signals were 
transmitted at different times, so that the term containing these clock errors describes 
the change in satellite clock error between transmit times. In theory, for the cesium 
clocks of the satellites this term may be considered as insignificant, although with the 
postulated future dithering of the signals this may well not be the case. In practice, the 
magnitude of this term depends on two factors: (a) the difference in transit times, and 
hence the length of the baseline between the two stations and their relative geometry 
with respect to the satellite; and (b) the level of synchronization between the two station 
time scales. 

Similarly, the equation of a between-receivers and -satellites, or double- 
difference carrier beat phase observable can be readily obtained from (4.23) for two 
stations (1 and 2) and two satellites (j and k) as 


<>j k 12 = <|> k i2- <t> j 1 2 (4.24a) 

• • 

= - (f/c) {[p k (T 2 ) - p k (T 1 )] + [5 t 2 P k (T 2 ) - Sx, p k (T i )]} 

- f [dt k (t k 2 ) - dt k (t k , )] + f [dT 2 (T 2 )-dTi(Ti)] 

-[ 6 k 2 

,atm ‘ 5 k l,atml + [N k 2 - N k n ] 


- { - (f/c) «pi(T 2 ) - pi(T 1 )] + [5x 2 pi(T 2 ) - St, pid,)]} 
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- f ldti(ti 2 ) - dtJ ( tJ i )] + f (dTgdgJ-dT^T,)] 
-[Sjg.atm-Sjl.atml + ^-NM} 


(4.24b) 


= * (f/c) {[p k (T 2 ) - p k (T i )] - tpi(T 2 ) - pi(T 1 )]} 


- (f/c) {[5 t 2 pk(T 2 ) - 5^ pkdO) - [5 x 2 pi(T 2 ) - St, pi(T,)]) 

- f «dtk(tk 2 ) - dtk(fk 1 )] - [dti(ti 2 ) - dtl(ti 1 )]} 

* {$ k 2,atm - 5 k i atm ] - [5j 2 atm - 8J'i ia tm]} 

+ [N k 2 - N^] - [Ni 2 - NJ-, ] . (4.24c) 


In this equation ail receiver clock errors have been removed, whereas the term 
containing the satellite clock errors has been reduced by a factor c/p =1 0 5 compared to 
that of the single-difference equation, thus effectively making this observable also 
insensitive to the satellite clock errors. 

Double differences are mathematically correlated because of the differencing 
process. This correlation depends on how the double differences are formulated from the 
corresponding single differences. For instance, if they are always formed between two 
successive satellites, the correlation coefficient between two so-formed double 
differences will be -0.5. Preferable to this scheme, is to account for this correlation 
(as it was done in the software of the present study) by considering the measurements at 
the same time epoch and to construct an elementary transformation 


Cqq = D Cgo D t 


(4.25) 


where Cgo is the covariance matrix of the individual single differences at that epoch 

assumed uncorrelated and with variance a 2 ; D is the double-difference operator 
transforming the single-difference observation space into the double-difference 
observation space, i.e. 
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-10 0 


0 ... 0 0 

0 ... o 0 


(4.26) 


D = 


1-1 0 

1 0 0 0 ... - 1 0 

1 0 0 0 ... 0 - 1 


leading to the covariance matrix associated with these double-difference observations 


C 

CD 


a 2 


21 11 ...11 

12 11 ...11 

11 11 ...2 1 

11 11 ... 1 2 


(4.27) 


or equivalently to the weight matrix 


P 

DO 


q~ 2 
n + 1 


n - 1 - 1 - 1 ... - 1 - 1 

- 1 n - 1 - 1 ...- 1 - 1 

- 1 - 1 - 1 - 1 ... n - 1 

- 1 - 1 - 1 - 1 ... - 1 n 


(4.28) 


which is of dimension n x n where n+1 is the number of satellites observed at that 
particular observing epoch. The weight matrix for an entire session has a block- 
diagonal form with each individual sub-block being of the form (4.28) and the size of 
different sub-blocks being generally variable depending on the number of satellites 
observed at that particular epoch corresponding to each sub-block. 
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SECTION 5 


STRATEGIES FOR ORBIT DETERMINATION 
AND BASELINE ESTIMATION 


The objectives ot future GPS-based geodetic systems have often been summarized 
as being twofold: 

- to provide accurate positions for the thousands of horizontal and vertical 
control points that comprise the current national geodetic databases in support of the 
uniform topographic mapping efforts and the growing need for better geodetic control in 
multipurpose cadastre systems; this would require both stations for remote (100’s of 
kilometers) positioning with decimeter accuracy, and stations for local (10’s of 
kilometers) positioning to a 2-3 cm accuracy level; 

- to measure vector baselines over distances from 100 to 4000 km both in 
length and orientation with an accuracy better than 0.01 ppm and to apply this 
capability to the study of the tectonic and crustal deformations. 

To fulfill these objectives, and in order to fully utilize available time in the field 
without an increase in the workforce, it is desirable to have highly compact stations 
which can be easily transportable, set up and operated by a minimum number of 
personnel (e.g. 1-3 persons) and by implementing automated techniques, be able to 
carry out such measurements over short observation periods (e.g. less than 4 hours 
including calibration assemblies). Any conceivable strategy to meet these needs must be 
three-tiered: the highest level being the orbital positions of the GPS satellites, the 
intermediate level being a set of continuously operating GPS receivers at several 
kilometers’ spacing, and the lowest level, the traditional monumented control. There are 
currently two independent efforts leading to the development and setting up of such a 
general system : 

- the development by the Jet Propulsion Laboratory for NASA (and independently 
by the U.S. Geological Survey (USGS), the National Geodetic Survey (NGS) and the 
University of Texas) of a general system based on the Fiducial Network concept; and 

- the development of an Active Control System (ACS) currently being 
implemented by the Canadian Geodetic Survey and the Geological Survey of Canada and in 
some slight variation in scope by the Texas Department of Highways and Public 
Transportation (Merrell, 1986). 
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5.1.1 The Fiducial Network System 


The Fiducial Network System (FNS) development, like the ACS (to be briefly 
described below), leads to a design that includes three major components: mobile GPS 
receivers, a network of well known fiducial stations and a central processing facility. 
The GPS data acquired simultaneously at the mobile and fiducial sites are brought to a 
central processing facility where the differential observables are reduced, calibrated 
and analyzed to jointly determine the orbits and the geodetic baselines in the framework 
of the few fiducial sites whose coordinates are accurately known from VLBI or SLR 
(Davidson et al., 1985). The need for cm-level baseline accuracies over distances of 
several thousand kilometers imposes a requirement for a regional network of fiducial 
sites whose relative positions are known to approximately to 1-3 cm (Dixon, Golombek 
and Thornton, 1985) which can only be provided by VLBI or SLR and improved (by some 
factor of 3) tropospheric calibration accuracies. The extent of these and other 
calibration problems is still the subject of investigations by JPL and several other 
groups. Obviously, for some applications this requirement may be relaxed somewhat to 
a few cm considering that the sensitivity of the results to fiducial baseline accuracy is 
nearly proportional to the length of the baseline to be measured. 

The same cm-level requirement in baseline accuracy from a few hours of 
observations translates into a requirement of a few cm of measurement noise over a few 
seconds including system noise, uncalibrated atmospheric effects (hence, simultaneous 
acquisition of the Li and L 2 carriers to assist in the ionospheric calibrations, and 

preferably the Li and L 2 code to assist the ambiguity resolution, is a critical support 
requirement), multipath, and other instrumental errors. Similarly, the requirement of 
simultaneous tracking of all satellites in view dictates the use of omnidirectional 
antennas or appropriately designed phase array antennas. 

The central processing facility for FNS should be suitably equipped to receive and 
process all data taken by the fiducial and mobile site receivers, determine precise orbits 
and post-fit GPS ephemerides self-consistent over extented periods at least at the level 
of 0.1 ppm or better. 

To satisfy the minimal regional geodetic requirements, at least three fiducial 
sites are required in order to provide substantial support for the determination of 
orthogonal (i.e. east-west, north-south) baselines. Ideally, the fiducial sites should be 
sufficiently close to the areas to be surveyed so that adequate simultaneous viewing 
between mobile and fiducial sites is assured, whereas the mobile sites can lie inside or 
outside the fiducial network (as long as they are within the effective area of the FNS 
sites) without suffering noticeable degradation in recoverable baseline accuracies. 
Theoretically, this range of mobile stations from FNS sites can be of the order of 
thousands of kilometers due to the high GPS orbits, however, in practical terms the 
optimal effective range will be somewhat limited by the correlation distance of 
atmospheric conditions affecting refraction, i.e. typically 200-300 km. 
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Fig. 5.1 - Schematic representation of the Canadian GPS-based 
Active Control System station operation. 
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5.1.2. The Rctiue Control System 


The Active Control System (ACS) under development by the Canadian Geodetic 
Survey (CGS) is based on the concept of a sparse network of control stations where 
continuous GPS tracking is carried out. A number of these Active Control Points (ACP) 
are to be permanently collocated at a number of widely spaced, fixed VLBI antennas thus 
ensuring the highest possible stability and long-term continuity of the ACS framework 
and consistent scale and orientation of the terrestrial networks positioned by GPS. With 
continuous tracking capabilities, the ACS will perform four major functions: (a) the 
data will be used along with similar data collected at stations like the FNS sites in other 
countries for the computation of precise orbits for distribution to various GPS users; 
(b) user control network connections based on ACP and user data sets combined. Because 
the coordinates of the ACS stations will be precisely determined (due to the continuous 
collection of data in these stations), they will provide precise control without the 
requirement for GPS users to occupy control stations; (c) broadcasting differential 
corrections for navigation purposes; and (d) as a means of verifying the integrity of the 
GPS system. Currently, the VLBI station at Yellowknife is the first such prototype ACP 
station in Canada in operation, which is also contributing data to the NGS/UoTexas 
Fiducial Network System between Austin, Texas, Westford, Mass, and Mojave, Ca. At 
least two new stations are planned to be brought "on-line" as soon as early next year as 
part of the ongoing developments. For the final ACS system a complement of as many as 
12-14 stations are envisaged (Delikaraoglou et al., 1986). 

The ACS system is being designed with 4 major features: automatic tracking 
stations, a tracking system design capable of collecting both code and carrier beat phase 
data, the capability to generate precise station position estimates with data taken over 
short periods of time, and the capability of the system to "self-check" its performance. 

In Figure 5.1 the major elements of the system are depicted in a functional form. 
An Active Control Point (ACP) unit would consist mainly of GPS receiving equipment 
(capable of tracking all visible satellites at the site) and a microcomputer housed in 
either existing facilities (e.g. provincial Mapping Agencies) or in unattended, 
environmentally-controlled shelters. Automatic weather sensors and likely a high- 
precision frequency standard will also be part of the ACP units. Currently, in the 
prototype stage of the development, TI-4100 receivers are used to evaluate the various 
functions of the ACP operations. However, the choice of GPS receivers to be eventually 
deployed in these sites is not planned to be restricted to one particular type. 

The planned operations of the ACPs require that all the functions of the system 
are microcomputer-controlled on the site with most functions such as satellite receiver 
control, tracking schedules, data acquisition tasks, data management, etc. be 
accomplished by software residing in the controlling microcomputer without 
interruption of the receiver tracking functions. Likewise, other housekeeping functions 
at the ACP site such as tracking status, error diagnostics, station operational status, and 
recording of weather data from the meteorological sensors will also be carried out by the 
ACP controller. 

The data collected at each ACP site will be appropriately edited, verified and 
suitably compressed before it is stored locally on disk or transferred via the chosen 
communication link to the Master Active Control Station (MACS). MACS is essentially 
the central control node for the ACS network and the control center of the ACS’s 
continuous operations and quality assurance. It is also the facility where computation of 
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GPS precise orbits would take place. Effective real-time data distribution on a regional 
and continental basis will depend on cost-effective mobile satellite communication 
services being available. This is expected early in the 1990's with the use of the new 
satellite communications system known as the Mobile Satellite (MSAT) which is 
currently entering its final stage of development and is expected to be fully deployed 
about the same time as GPS becomes fully operational and the full ACS network is in 
place (Delikaraoglou, 1986). 


5.2 Modes of Analysis for Orbit Determination and Baseline 
Estimation 


The problem of accurately modelling the forces acting on the GPS satellites has 
been discussed to some extent in Section 3 where it was indicated that forces which 
contribute to the satellite motion comprise: 

- gravitational attractions of the earth; 

- gravitational attractions of the sun, moon, and the planets (i.e. the third-body 
effects); 

- the variable part of the earth’s gravitational fields arising from the tidal and 
other deformations of the solid earth and the oceans; 

- solar radiation pressure (both direct and albedo effects including the so called 
Y-bias effect (Fliegel et al., 1985)); 

- atmospheric drag forces; and 

- magnetic forces and other smaller effects. 

The central problem in any orbit determination scheme is to accommodate a 
relatively sophisticated model for the above force fields with the degree of sophistication 
being, loosely stated, a function of the length of the arcs under consideration. Whereas in 
theory, as already seen in Section 3, the solution of this determination problem involves 
a rather complex integration process of the satellite’s equations of motion, in practice, 
it ultimately reduces to an orbit improvement problem using observations which are 
themselves functions of many satellite positions, ground positions and other nuisance 
parameters such as satellite and receiver clock errors, atmospheric delay parameters, 
etc. In the case of the GPS satellites and the carrier beat phase measurements, this 
usually involves a simultaneous adjustment of: 

- satellite epoch states; 

- station and satellite clock parameters; 

- selected station locations; 

- zenith tropospheric delays; 

- satellite solar radiation pressure coefficient(s); 

- empirical accelerations; and 

- carrier phase ambiguities. 

Over the past 2 years high-precision orbit determination results for the GPS 
satellites have been reported (e.g. Abbot et al., 1986; Beutler et al., 1986; Williams, 
1986) from which the quality of baseline results have improved by nearly two orders of 

magnitude over that period; i.e., from several parts in 10 6 in 1984 to between 2 parts 

in 10 7 and 5 parts in 10 8 today for the longer baselines. In achieving these results, two 
methods of orbital analysis have been used almost exclusively: 
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- the "Fiducial Network" mode based on the approach just outlined in which the GPS 
orbits and determined baselines are defined in the framework of a few fiducial stations 
whose coordinates are accurately known from VLBI or SLR; and 

- the "Free Network " mode approach (e.g. Beutler, 1985; Delikaraoglou, 1985) in 
which a combination of station and satellite coordinates are allowed to vary 
simultaneously in an adjustment by way of imposing relatively strong a priori 
constraints on the orbits and letting the terrestrial network adjust free. 

Fundamental to each of these analysis techniques are also considerations 
regarding the length of the satellite arcs vis-6-vis the network configuration. Single- 
pass arcs (about 1/3 satellite revolution or 4-5 hours duration) require less elaborate 
force models (cf. Buffett, 1985) and hence, fewer parameters to solve for, whereas 
longer multi-day arcs (1 day to 1 week duration or more) are preferable in order to 
define a more stable reference frame, but they require additional care to account for any 
unmodelled (or improperly modelled) contributions of the forces acting on the satellites. 
The latter are particularly important for long arcs since, if unaccounted for, they can 
adversely affect the results in a systematic way thus becoming (along with data noise) 
the limiting error source. Tracking networks of different scales and geometries allow us 
to differentiate between orbital and other distance-dependent error sources in GPS such 
as unmodelled atmospheric delays. To date for instance, the U.S. Department of Defense 
has relied on its global network of the GPS monitor stations and the use of the 
pseudorange and integrated Doppler observables to produce predicted (broadcast) orbits 
at the tens-of-meters accuracy level. However, for the more stringent orbital 
requirements of precise geodetic and geodynamic studies with GPS, the optimal network 
configuration for the determination of GPS orbits and station coordinates is not well 
established yet. Ideally, dense regional and local networks in the area of interest at the 
scales (i) 100-500 km, (ii) 500-2000 km, and (iii) 2000-4000 km would provide 
the test ground for isolating possible existence of systematic errors specific to the 
space/terrestrial conditions at the time of the measurements and thus would allow for a 
comprehensive study of various orbital analysis techniques, particularly the ones that 
use increasingly shorter satellite arcs. Results of various comparison tests carried out 
as part of this study to assess each strategy’s practical capability to reduce particular 
systematic errors and improve the orbit and baseline determinations will be reported 
later on in this section. Before doing so, however, we will discuss another outstanding 
item tacitly assumed as part of each of these modes of analysis: the effect of the a priori 
covariance information on the so-determined orbits and baselines. 


5.4 fl priori Information In the GPS Satellite Netmorks 


In any adjustment of differential GPS observations the various unknown 
parameters can conveniently be grouped into three categories: 

- station parameters or baseline components; 

- satellite epoch state vectors; 

- other dynamical satellite parameters like solar radiation pressure, 
instrumental-type parameters like clock errors and other non-geometric 
parameters like atmospheric errors. 
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In practice, provided that enough differential observations of adequate geometry 
of adequate geometry have been collected, a least-squares adjustment can be carried out 
in which the primary interest is in the coordinates of the terrestrial points and hence, a 
solution is often sought from a partitioned system of normal equations of the form 




(5.1) 


where 


= [ A c | A b ] 

first partitioning 

(5.2a) 

= [Ag A s 1 A b ] 

second partitioning 

(5.2b) 

= ( Af A m | A s | A b ] 

third partitioning 

(5.2c) 


is the partitioning of the design matrix A of the network; 

X = [ X c X b ] T first partitioning (5.3a) 
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= [XgX s I X b ]T 
-[X, x m I X s I X b ]T 


second partitioning 


(5.3b) 


third partitioning (5.3c) 


is the corresponding partitioning of the unknown vector into a sub-vector X c of 
coordinate unknowns (both ground X g and satellite X s ) and a sub-vector X b of nuisance 
parameters (non-geometric, instrumental, etc); X g is further partitioned in a set of 
coordinate unknowns Xf corresponding to fixed or constrained (i.e. fiducial) stations and 
X m corresponds to mobile or new stations to be surveyed; Njj are, respectively, the 

partitioned parts of the normal equations pertaining to the partition of X; Pj = C'^Xj) 

are similarly weight matrices corresponding to the a priori information (i.e. C'^Xq j) 
the a priori covariance matrices) available on the parameters X 0 assuming that there 
is no correlation between the different sets of parameters. 

The correspondence of this partitioning scheme to the aforementioned two main 
modes of analysis is now readily apparent. The "fiducial network" mode corresponds to 
the case whereby 0 < P s < °°, Pf -» °° and 0 < P m < «. Respectively, the "free 
network" mode case corresponds to 0 < P s < °o, 0 < Pf < °° and 0 < P m < <» whereby 

furthermore, the first two sets of conditions of a priori information are necessary since 
they generally reduce the datum defect inherent in the free network adjustment 
(Delikaraoglou, 1985). Generally there are a number of practical consequences arising 
from the a priori information on either the satellite parameters or the station 
coordinates of some of the stations and the effect which they may have on the variances of 
the coordinates of the rest of the determined stations and the baseline components 
between such stations. These can be readily understood by examination of the more 
general form of the normal equations (5.1). 

It can be shown that starting from the second partitioning of the parameter 
vector X, and after successive elimination of X b and then X s and some lengthy 

derivations, the post-fit covariance matrix of the station coordinates X g including a 
priori information is given by 


C(Xg ) = 


( N + P ) - N ( N + P I’’ N T 

ff f fs SS s 


fs 


N . N ( N + P ) * 1 N T 
fm fs ss s ms 


N T . N ( N + P ) 1 N T 
fm ms ss s f s 


N . N ( N + P ) 1 N T 
mm ms ss s ms 
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where the submatrices Qjj are given by 

Qi 1 = [ B i 1 - B i2 ( b 22) _1 ( b 12) T l* 1 

- I K Nff + Pf ) * Nfs ( N ss+ P s )-1 (N fs ) T ] 

- (( Nfm - Nfs ( N ss + P s ) 1 (^ms ) T ) 

X I Nmm ' N m s ( N ss + P s )'^ (N ms )^ J" 1 

* {(N fm ) T - N ms ( N ss + P s )-1 (N, S ) T } ]-1 (5.5a) 

Q 2 2 = [ (B 22)' 1 * ( B 22) _1 ( B 1 2) T Ql 1 B 12( B 22) _1 1 
= [B 22 -(B 12 ) T (B 11 )-1 b 12 j-i 

= [[( ^mm ) ‘ N ms ( N ss + P s )’1 (N ms )T ] 

- (( N fm ) T - N ms ( N ss + P s )-1 (N f8 )T } 

X [ (N ff + P f ) - N fs ( N ss + P s )-1 (N f8 )T ]-1 

x { N fm - N fs ( N ss + P s )-1 (N ms )T } ]-1 (5.5b) 

Qf 2 = - Qfl B 12 (B 22 ) _1 

= [B 22 -(B 12 )T( Bii )-1 b 12 ]-i 
= -Qff [ N fm - N fs ( N ss+ P s )-1 (N ms ) T J 

X I (N mm - N ms ( N ss + P s ) _1 (N ms )T ]‘1 (5.5c) 

0 2 1 = (Ql2 ) T • (5.5d) 
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The covariance matrix for the baseline components between stations in the f- and m-sets 
respectively are readily obtained by a linear transformation: 



Q 

Q 

11 

12 

T 

Q 

Q 

12 

22 


” Qi i + Q 22 • Qi 2 ' (Qi 2 ) T 


- Qll + Q22 + Q11 B12 (B 22)' 1 
+ (B22)- 1 (Bi 2 ) t Q11- 


(5.6a) 


(5.6b) 


(5.6c) 


The effect of a priori information can now be readily seen for the following cases: 

(a) the satellite parameters are constrained and there is no knowledge on the coordinates 
of the f-set of stations (P s -4 00 , p f 0 ), e.g. in a case of reducing GPS data at 

completely unknown stations with orbits assumed completely known and held fixed in the 
adjustment. Under these assumptions equations (5.4) through (5.6) reduce to 

C(X,) - [ N„ - N lm (N mm )-1(N, m )T 1-1 (5.7a) 

C(X m )-[N mm - (N fm )T( Nf ,)-1(N, m ) ]-1 (5.7b) 


C(B,, m ) . C(X f ) + C(X m ) ♦ C(Xf) (N, m )T (N mm )-' 


+ <N mm )-’<N (m ) T C<X,); (5.7c) 

(b) the satellite parameters are partially known and the coordinates of the f-set of 
stations constrained (0 < P s < «>, Pf -» °°), e.g. in the case of reducing GPS data from 

fiducial stations for the determination of orbits along with the coordinates of unknown 
stations (m-set of stations) 


C(X f ) = 0 (5.8a) 

C(X m )-[N mm -N ms (N sst P s )-'(N ms )T]-1 (5.8b) 
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C(B f , m ) = [ N mm 


- N ms (N ss + Ps^Nr^s) 1 1' 1 = C(X m ) ; (5.8c) 


(c) satellite parameters and the coordinates of the f-set of stations constrained 
(P s -» oo, Pf oo), e.g. the case of reducing GPS data between fiducial and unknown 

stations and holding the orbits (as determined from the fiducial stations, for instance, in 
a previous run) fixed during the determination of the coordinates of the unknown 
stations (m-set of stations) 


C(Xf) - 0 

(5.9a) 

C(X m l - 

(5.9b) 

C(Bf. m ) - (Nmm)' 1 - C<X m ) ; 

(5.9c) 


(d) satellite parameters and the coordinates of the f-set of stations partially known 
(0 < P s < oo, 0 < Pf < oo), e.g. the case of reducing GPS data in the free network mode by 

imposing relatively strong constraints on the orbits and some stations (i.e. the f-set of 
stations) and letting the rest of the network (i.e. the m-set of stations) free to adjust 


C(Xf) 

» equation (5.5a) 

C(X m ) 

= equation (5.5b) 

C(B,,m) 

- equation (5.6). 


A careful examination of these covariance expressions leads to some useful 
insight as to how the results may vary under different circumstances of usage of a priori 
information. For instance, from (a) and (c) it is readily seen that when the satellite 
orbits are assumed known, the uncertainty of the coordinates of new surveyed stations 
increases as the uncertainty of the fiducial stations increases. A similar increase is also 
to be expected in the new stations when the coordinates of the fiducial stations are 
constrained and the uncertainty in the satellite orbits increases; i.e., case (c) vis-a-vis 
case (b) and a further increase is to be expected when both the uncertainty in the 
satellite orbits and the coordinates of the f-set of stations increases; i.e., case (b) vis- 
a-vis case (d). The last observation indicates a significant weakness, and as it turns out 
from a practical viewpoint, a problematic area of the free network mode of analysis vis- 
a-vis the fiducial network mode: unless appropriate, realistic orbital constraints are 
chosen for the former, variable results in the solution are to be expected. This fact alone 
suggests that the free network mode should be an unlikely candidate to use for precise 
geodynamic applications where the utmost certainty in the results is a must. 

For the baseline components, the examination of the covariance expressions 
(5.8c) and (5.9c) indicates that given that the coordinates of fiducial stations are often 
assumed known, increasing the uncertainty in the satellite orbits leads to an increased 
uncertainty in the components of the baselines connecting fiducial and non-fiducial 
stations; also if the satellite orbits are assumed known, increasing the uncertainty of the 
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fiducial stations also increases the uncertainty in the components of the baselines 
between fiducial and non-fiducial stations (cf. Eq. (5.7c) vs. Eq. (5.9c)). 


5.4 EHperimental Results 


The methodology described thus far has been implemented in software and tested 
with various GPS data sets. In this brief overview, we shall report on representative 
tests relating to different aspects of orbit determination and the results obtained with 
sample data available from two such GPS campaigns designed to deploy a large number 
and distribution of receivers (Fig. 5.2). The first campaign took place in the Spring of 
1985 and involved TI-4100 receivers (Henson et al., 1985) deployed at some ten sites 
in Southern California, JPL Series-X receivers (Crow et al., 1984) operating at 
Mojave and Owens Valley in California, and prototype Macrometer*ll receivers (Ladd 
and Councelman, 1985) operating at the three POLARIS sites at Haystack, MA, 
Richmond, FL and Ft. Davis, TX. Six days of data (days 90-95) from this campaign 



Fig. 5.2 - Site locations for the Spring '85 and Summer '86 GPS Experiments 


were used in our tests. The second campaign took place in the Summer of 1986 for the 
main purpose of establishing first epoch measurements in Iceland, but it also involved 
observations from GPS receivers providing simultaneous coverage of the northern 
hemisphere in Canada, USA, Sweden and Denmark. For the purpose of this report, 
results with only a 2-day portion of the data (days 193-194) collected in Northwestern 
Canada during this campaign will be reported. 
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Most of the analyses reported herein were carried out using a modified version of 
NASA’s GEODYN precision orbit determination program and software developed for the 
preprocessing of the raw GPS data (including formatting, editing, calibration and cycle- 
slip detection). GEODYN provided integrated trajectories in which the earth’s gravity 
field model was expressed by spherical harmonics (up to degree and order 8 for the case 
of the GPS satellites in accordance with the discussion in Section 3); point mass 
attractions from the sun and moon and a direct solar radiation pressure model (i.e., a 
constant acceleration along the sun-satellite vector) were included, plus the model of 
empirical accelerations developed by Colombo (1986), i.e. Eq.(3.36). 

The observations used to compute the orbits for the Spring ’85 Experiment were 
from Macrometer data collected at the three POLARIS sites. A copy of these observations 
was transferred already "cleaned" from cycle slips with the compliments of the MIT GPS 
group. Data from the TI-4100 receivers at Owens Valley, Mojave, Ft. Davis and Hat 
Creek were pre-processed at GSFC mostly with software developed in the course of this 
study. Similarly, TI-4100 data from the Summer '86 GPS campaign were transferred 
from CGS and pre-processed at GSFC for editing, phase connection, data compression and 
atmospheric calibrations. 

The observables used in our analyses were basically "between-stations" single 
differences. Dual-frequency and L 2 observations were used throughout. Since this 

type of observable is affected by receiver clock errors, the clock behaviour was 
modelled as "white noise"; i.e., one common clock bias was solved for per measurement 
epoch for stations with no external frequency standards; this is effectively equivalent to 
double differencing but it requires many more bias parameters-an actual double 
differencing scheme is currently implemented in GEODYN which hopefully will eliminate 
the constraints imposed by this equivalent "white noise" model in the future. 

A number of tests were carried out to investigate the differences between the 
fiducial concept and the free network approach in both the orbit determination and 
relative positioning with GPS vis-&-vis the impact on these determinations of the length 
of the satellite arcs and the variations in the tracking network configurations. For this 
purpose identical data sets were processed in each of these modes and the results were 
compared. 


Spring ’85 Experiment tests 


For the fiducial network strategy, typically, the coordinates of the stations at 
Haystack, Richmond and Ft. Davis were held fixed at their VLBI values, whereas the 
coordinates at Owens Valley were adjusted along with the satellite orbits, clock biases, 
zenith tropospheric delays, etc. This configuration provided a well-distributed set of 
points to test the orbit quality and baseline repeatabilities over distances of 1000- 
4000 km. Descriptions of various results obtained from the analysis of this data set 
were also reported by Pavlis et al (1987) and Colombo and Zelensky (1987). 

For the free network approach, Haystack was held fixed and the coordinates of all 
the other stations were estimated along with the satellite orbits. As already mentioned, 
this approach requires that some a priori (variance/covariance) information about the 
orbits be introduced into the normal equations of the adjustment process. In the results 
to be reported below, a priori rms errors for the state vectors were chosen at the rather 
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Fig. 5.3 - Orbit precision with 6-day multi-arc solutions 
using TI-4100 data from the March '85 Experiment 


extreme values of 20 m in position and 0.2 m/sec in velocity. Several tests in the large 
domain 2-20 m and 0.2-2 m/sec were also carried out with variable results depending 
on the imposed constraints. As already pointed out, this is currently one of the most 
problematic areas of this approach, that is to say, to come up with realistic estimates of 
the imposed constraints on the orbits. This also verified the practical problem hinted in 
the examination of the covariance expressions in the previous section and confirmed the 
adverse effect in the solution results under different circumstances of a priori 
information. 

Fig. 5.3 shows plots of formal orbit errors from a 6-day solution for the GPS 
satellites 6, 8, 9 and 1 1 indicating the level of precision of the computed orbits. GPS 1 1 
was a well-tracked satellite (i.e. with good temporal coverage) during this period, 
whereas the other satellites were more sparsely tracked each day. As the figure 
indicates, in the fiducial network approach orbit precision of better than 1 m in all 
position components is achievable from multi-day arcs for a well-tracked satellite, and 
2-3 m for sparsely tracked ones. By comparison, for the free network approach, the 
formal errors are consistently higher for GPS 11 in spite of its good tracking coverage, 
whereas they are quite comparable for the other three satellites (somewhat worse 
results were obtained in another test whereby Ft. Davis was held fixed instead of 
Haystack in order to test the sensitivity of this approach to the choice of the fixed 
station). 

Since in general, formal errors would tend to underestimate orbit accuracy, a 
comparison between solutions of two overlapping arcs of different length were made in 
which the recovery of the position and orientation of Owens Valley as determined with 
respect to the three POLARIS sites was compared. The discrepancy vector obtained in 
each case was resolved for an easier comparison into a coordinate system whose axes are 
aligned 

- along the true baseline (in this case the VLBI baseline) from the fixed to the 
adjusted station; 

- along the axis in the direction of increasing azimuth for the baseline; and 

- along the axis in the direction of increasing elevation for the baseline (which 
is almost identical to the local geodetic height component). 
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Fig. 5.4 - GPS vs. VLBI at Owens (non-fiducial station) illustrating 
the effect of variation of the arc length 


Fig. 5.4 shows the results of this comparison for a 2- (days 90-91) vs. a 6-day 
arc (days 90-95) in each of the fiducial and free network modes respectively. In the 
fiducial network case, the azimuth discrepancies of the three baselines Haystack/Owens, 
Richmond/Owens, and Ft. Davis/Owens, as well as the discrepancies in their lengths 
(3929, 3741 and 1508 km respectively) are all within 10 cm for the 2-day arc, 
whereas the elevation discrepancies are all at the 25-cm level with respect to the VLBI 
values at Owens. The corresponding azimuth discrepancies in the free network case are 
also within 10 cm, whereas the baseline length and elevation discrepancies show a 
reverse tendency compared to the previous case. That is, elevation discrepancies were 
reduced to the 15-cm level, whereas baseline discrepancies up to 20 cm were noted. As 
may have been expected, there was a very noticeable improvement in all cases when the 
6-day arc results were compared. The fiducial network results were for all baselines 
within 5 cm in azimuth and length, and 10 cm in elevation which corresponds to better 

than 1 part in 10 8 relative to the length of these baselines. Considering that the 
accuracy of the GPS orbits and relative positioning accuracies are related by the rule of 
thumb 


(dr/r) = (dB/B) 


where B is the baseline length, r is the average range to the satellite (about 26,000 
km), dB is the baseline error and dr is the orbit error, a 10-cm length discrepancy (as 
noted in the fiducial network case) over the 3000-km average length of these baselines 
would infer an orbit error of 0.8 m. Similarly a discrepancy of 20 cm (as noted in the 
free network case) would correspond to 1 .7 m orbit error. These values may give a more 
representative measure of the orbit accuracy than would formal errors and would 
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indicate what one could expect from multi-day arcs for well-tracked satellites in each 
case. 


In order to assess further the quality of the so-computed GPS orbits, a set of 
daily solutions was carried out to test baseline repeatability. Using TI-4100 data, the 
station locations at Mojave and Hat Creek with respect to Owens Valley were 
independently computed each day. The orbits used in these solutions were those computed 
from the 6-day arcs using Macrometer data from the three POLARIS sites and Owens 
Valley and they were held fixed in these tests, so these comparisons provided another test 
not only of the accuracy of the orbits but also of their robustness, and of the models used 
to provide the integrated trajectories. The noted differences from the VLBI values are 
shown in Fig. 5.5a and 5.5b. Four daily solutions and their respective means and 
associated standard deviations are shown for the east, north and vertical directions and 
the corresponding length of each baseline is also shown. In the case of Mojave, the 
location of the SLR monument is also indicated as a further measure of comparison. In 
the case of the Owens-Mojave (245 km) baseline the standard errors of the mean in all 
components except the height were better than 5 cm which corresponds to a precision 
better than 2 parts in 10 7 . 
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Fig. 5.5 - Day-to-day repeatability (a) above, Owens-Mojave (245 km) 
(b) below, Owens-Hat Creek (484 km); 6-day orbits (days 90-95) 
kept fixed in daily solutions. 
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The higher scatter observed in the respective height differences corresponds to a 
precision of 5 parts in 10 7 which is probably due to the troposphere. It would appear 
that for the shorter baselines, tropospheric delay fluctuations may be affecting 
repeatability more than is the case with the continental baselines where orbital 
inaccuracies may be the more dominant source of error. From the North/East plot, it is 
interesting to note that there seems to be a great deal of correlation between the noted 
differences and the direction to Owens which may be indicative of the geometric weakness 
of the daily observability windows for this baseline. Similar results were obtained for 
the daily solutions carried out for the Owens-Hat Creek (484 km) baseline. Again all 
components, including the height, showed a scatter indicating a precision at the level of 2 

parts in 10 7 - A similar correlation of the discrepancies in the north and east 
components and the direction to Owens (note that Hat Creek, Owens and Mojave are 
nearly along the same line) is also noticeable. 


Summer *86 Northwestern Canada Campaign tests 


The tests performed thus far were with a valuable body of data available from 
sites with well-known coordinates from VLBI and from several days of observations. 
There are situations however, when none of these conditions are fulfilled. In such a 
situation the free network approach can be proven extremely useful. This was the 
situation, for instance, with the 2 days of data collected at the regional network in 
Northwestern Canada in the Summer of 1986. Although Yellowknife, Whitehorse and 
Penticton are VLBI sites in Canada, the connections of the GPS markers to their VLBI 
counterparts were not available which rendered their use as fiducial locations useless 
for the purposes of our tests. It was possible however to conduct further tests using the 
free network approach. Here we just give an outline-without presenting complete 
results (which are to be reported in a subsequent detailed report)-of the evaluation of 
this data set using orbits computed with initial conditions derived from the broadcast 
ephemerides. 

On days 193-194 (July 11-12) four TI-4100 receivers were operated by CGS 
at the above sites and the Pacific Geoscience Center (Pat Bay, Sidney, B.C.). Satellites 3, 
6, 8, 9, 11, 12 and 13 were observed on both days for 3-5 hours each night, in 
addition to these satellite passes these satellites were observed 12 hours earlier, but 
this data was not incorporated in these tests. In summary, our processing consisted of 
the following steps: 

- We used the broadcast ephemerides to obtain the initial conditions needed in the 
CGS software (in which similar models to the ones used in GEODYN have been 
implemented) to integrate the orbits: 

- In these integrated trajectories a force model similar to the one used in the 
previous tests was used, except for the solar radiation model which included only one 
radiation pressure coefficient per day which was considered adequate for the short arcs 
in these tests; 

- The Li and L 2 observations were combined to form the ionospheric-free linear 
combination: 

- Clock errors were similarly treated as "white noise"; and 
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• The observations were processed individually for each day using the following 
options: 

( 1 ) the coordinates of Yellowknife were constrained to their VLBI values with a 

priori variances of 1 m 2 per coordinate; 

(2) the coordinates of Whitehorse, Penticton and Pat Bay were treated as 
unknowns; and 
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Fig. 5.6 - Baseline repeatability using the "free" network approach 
with single pass, short-arcs for satellites 6, 8, 9, 1 1 . 


(3) satellite state vectors for each satellite arc were computed using 

constraints of 20 m in the Keplerian elements of each arc expect for the 
argument of the perigee value which were held fixed. 

It should be noted in this connection that in contrast to the multi-day arcs, the 
geometry is very limited in such arcs of a few hours’ duration. In this particular 
network-satellite configuration, for instance, there was considerable lack of data from 
the northwest direction which certainly compromised viewing geometry. Fig. 5.6 shows 
the baseline repeatability results using these single-pass arcs. The rms scatter of the 
two daily solutions of the components of the three baseline vectors (Yellowknife- 
Penticton/1 1 05 km, Yellowknife-Whitehorse/1495 km, and Yellowknife-Pat Bay/ 
1632 km ) about their weighted means were about 5 cm on the average and less than 15 
cm in the extreme cases noted in the baseline between Yellowknife and Penticton. 

Nevertheless these values are of the order of 5 parts in 1 0 8 which would suggest an 
orbit error of 1 .3 for these short arc orbits. 
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SECTION 6 


SUMMARY AND DIRECTIONS FOR FUTURE LUORK 


We have demonstrated, albeit with a limited data set, that dual-frequency GPS 
observations and interferometric-type analysis techniques make possible the modelling 
of the GPS orbits for several days with an accuracy of a few meters. The use of VLBI or 
SLR sites as fiducial stations together with refinements in the orbit determination 
procedures can and have greatly reduced the systematic errors in the GPS satellite 
orbits used to compute the positions of non-fiducial locations. In general, repeatability 

and comparison with VLBI of so-determined locations are of the order of 2 parts in 10 7 

to 5 parts in 10 8 for baseline lengths < 2000 km. The free network approach offers a 
powerful alternative technique capable of achieving nearly equivalent results for less 
than 2-day or shorter arcs (with somewhat worse results for longer arcs) particularly 
under conditions where one has no access to precise orbits or a network of VLBI/SLR 
fiducial stations in the vicinity of the survey. However, the requirement of having to 
define appropriate constraints needed to be imposed in the estimation process in this 
latter case is a drawback of this approach worth bearing in mind. 

Multi-day arcs (3-6 days) generally provide better orbits and baselines than 
single-day passes or up to 2-day arcs. In achieving these levels of repeatability and 
accuracy several refinements of the orbit error models are necessary. To that extent, 
Colombo’s (1986) general force model based on the likely resonant character of the GPS 
orbit errors seems capable of removing the main unmodelled (or incorrectly modelled) 
forces like solar radiation pressure, Y-bias, etc. thus confirming in practice the 
theoretical considerations that lead to its development. 

Considering the results of this analysis and similar results reported by several 
other groups, it seems that achieving high-quality orbits for the GPS satellites is, to 
date, a lesser problem than other error sources. From the results obtained and the 
experience gained so far it is clear that of the limiting factors to further orbit and 
baseline accuracy improvement, fiducial station coordinate errors, better handling of 
propagation medium errors (particularly of the wet troposphere), and further 
modelling of the orbital dynamics for the GPS satellites are areas to concentrate on for 
future work. All these errors can be better understood and reduced via tracking 
networks of various scales and geometries which would allow us to differentiate between 
the various error sources. The obvious question now is what levels of repeatability can 
be achieved with GPS by repeated occupations of the same sites at frequent intervals; the 
only reason for not having an answer yet is that not enough repeat experiments (other 
than the three occupations of some California sites during the Spring and Fall ’85 and 
again during the Summer ’86 campaigns) have been carried out so far. To date, several 
national and international campaigns (many of them NASA sponsored) have taken place 
and more are being planned for the future. Yet GPS work seems set to continue at many 
centers, in a largely decentralized way unless the data collected during these campaigns 
are shared for the purpose of carrying out such studies that hopefully will offer a unique 
opportunity to study these aspects further. 
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In retrospect, this report which is mainly a synthesis of problems, assumptions, 
methods and recent advances in the studies towards the establishment of a GPS-based 
system for geodesy and geodynamics marks only one phase in the continuing effort for the 
development of such a system. At the outset of this project, the physical phenomena 
relevant to the GPS carrier phase observable which does not require access to the 
(likely to be restricted in the future) coded information in the satellite signals, and 
general orbit force models were reviewed. Models for these phenomena were described 
for implementation in software, a large part of which was developed and tested in-house 
at GSFC. This activity, however, is by no means complete. Several enhancements and 
improvements are still possible especially if these models (and the software modules 
now available at GSFC) are to be eventually implemented in packages like NASA’s 
GEODYN II. Some of these include the full implementation of the double-difference 
module; the 6-parameter solar radiation general force model that does not require 
specific physical information on the satellites (e.g., panel shape, dimensions, material, 
etc.) and yet as demonstrated here is as good or better than those models that rely on the 
restricted information; adjustment strategies for simultaneous orbit and station position 
estimation; representation of refraction and clock parameters with time-varying 
polynomials or continuous stochastic processes of given correlation functions, and 
ability to write out residuals and selected partials for further processing outside 
GEODYN. Finally, it should be mentioned that although this study and the results 
reported here have relied on the use of the carrier beat phase observable, there might be 
good reason to experiment also with the synergistic mix of the carrier phase with the P- 
code pseudorange observable. The activity at JPL and recent results presented at the 
October 1987 Crustal Dynamics Meeting at NASA’s Goddard Space Flight Center have 
given positive indications that further improvements in precision can be achieved by the 
combination of the two. In any event, the sole fact remains that GPS results already 
agree impressively well with collocated VLBI results. To some, including the author, it 
now seems reasonable to expect within the next few years that more evidence will show 
GPS to be as powerful and reliable a tool as mobile VLBI and SLR are today, but largely 
more economical. Promises of pocket-sized receivers and wallet-sized prices in the 
not-too-distant future may well prove that GPS is indeed out of this worldl 
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